Chapter 5
Continuous random variables

Lectures 12–14 · measuring what happens: densities, the Gaussian curve, waiting times and the other shapes of chance

for measurements, single values have probability zero; intervals do not

Chapter 4 was about counting: how many passengers, how many fakes, how many calls. This chapter is about measuring: how tall, how long, how heavy, how concentrated. A person's height, the time until the next earthquake, the alcohol concentration in a blood sample, the age of a bone found at a crime scene. These quantities don't jump from one value to the next; they can take any value in an interval, and that changes the mathematics in a small but important way.

Here is the puzzle that drives the chapter. What is the probability that a randomly chosen man is exactly 70 inches tall? Not 70.1, not 69.99, but 70.000000... exactly. There are infinitely many possible heights, so the answer has to be zero. Yet every man has some height. How can every individual height have probability zero, when the probabilities must add up to 1? The answer is to stop asking about single values and start asking about intervals, and to measure probability as an area under a curve.

Here is the route through the chapter.

Sections 5.1 and 5.2 replace the PMF with a density curve, and redo expected value and variance with integrals instead of sums.
Section 5.3 meets the most famous curve in science, the Gaussian distribution, using Francis Galton's measurements of 481 Victorian men.
Section 5.4 is about waiting times: the exponential and gamma distributions, from mining disasters to radiocarbon dating.
Sections 5.5 and 5.6 meet the beta, log-Gaussian and other distributions, and Section 5.7 puts them side by side.

The calculus you need is in Section 0.4: areas under curves are integrals.

Densities serve the two goals of the course just as PMFs did. For prediction, a density fitted to past measurements tells us how often to expect future values in any range: how many sons will be over six feet, how long until the next mine disaster, how often an honest driver will be flagged by a breathalyser. For inference, we reverse the question: what does a measurement tell us about the hidden quantity behind it, such as a driver's true alcohol level, the age of a body, or a player's true shooting percentage? Along the way, watch for a third question that underlies both: has the mechanism that produced the data stayed the same? A model fitted to one era (the coal mines before 1890) can predict badly in another.

Think About It: Before You Start

A breathalyser reads 0.081 g/dL, just over the legal driving limit of 0.08. The device's readings are known to vary by a few thousandths around the true value. Should the driver be convicted? What would you need to know about the device to decide? Hold on to your answer: Example 5.7 works it out.

You'd want to know how big the measurement error typically is, and how likely a reading of 0.081 or more would be if the driver's true level were just below 0.08. If readings routinely vary by 0.005 either way, a reading of 0.081 is entirely consistent with a true level under the limit. Many jurisdictions deal with this by requiring two readings, or by subtracting a margin for error, and the size of that margin comes straight from the Gaussian distribution of measurement errors in Section 5.3.

Lecture 12 — Densities, expected value and the uniform

5.1 From histograms to densities

Galton's heights

In the 1880s the English scientist Francis Galton (a half-cousin of Charles Darwin) collected the heights of nearly 1,000 adult children and their parents, to study how height is inherited. His data still survive, and they're a good place to see how a continuous distribution behaves.

Example 5.1: Galton's Sons ★ lecture example

Galton's data include the heights, in inches, of 481 adult sons. (a) Draw a density histogram of their heights, with narrower and narrower bins. (b) What fraction were between 5 ft 7 in (67 inches) and 6 ft (72 inches)? (c) What does the histogram look like as the bins shrink?

The data are in the R package HistData as GaltonFamilies. To avoid installing anything, we read the same data straight from an online copy:

galton <- read.csv("https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/HistData/GaltonFamilies.csv")
sons <- galton$childHeight[galton$gender == "male"]
length(sons)
summary(sons)
mean(sons >= 67 & sons < 72)     # fraction between 67 and 72 inches

output

[1] 481
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  60.00   67.50   69.20   69.23   71.00   79.00 
[1] 0.6569647

(a) A density histogram (Section 2.1) scales each bar so its area is the fraction of sons in that bin. Here are three bin widths:

h <- rbind(data.frame(width = "2 inch bins", height = sons),
           data.frame(width = "1 inch bins", height = sons),
           data.frame(width = "half-inch bins", height = sons))
h$width <- factor(h$width, levels = c("2 inch bins", "1 inch bins", "half-inch bins"))
ggplot(h, aes(height)) +
  geom_histogram(data = subset(h, width == "2 inch bins"), aes(y = after_stat(density)), binwidth = 2, boundary = 0, colour = "white") +
  geom_histogram(data = subset(h, width == "1 inch bins"), aes(y = after_stat(density)), binwidth = 1, boundary = 0, colour = "white") +
  geom_histogram(data = subset(h, width == "half-inch bins"), aes(y = after_stat(density)), binwidth = 0.5, boundary = 0, colour = "white") +
  facet_wrap(~ width) +
  labs(x = "height of son (inches)", y = "density")

(b) About 65.7% of the sons were at least 67 and under 72 inches tall.

(c) Whatever the bin width, the outline has the same bell shape, peaking near 69 inches. With more and more data and narrower and narrower bins, the outline would settle down to a smooth curve. That curve is the density of the height distribution. The widget below does exactly this with simulated heights like Galton's: watch the bars smooth out into a curve, and watch the probability of a single bin shrink as the bin gets narrower.

A detail worth noticing. The half-inch histogram has spikes at whole inches. Heights were measured to the nearest quarter or half inch and many were written down as whole numbers, a common habit called heaping. That's why we asked about "at least 67 and under 72": an interval that cuts through a heap of recorded 72s would give a misleading answer. Measurement is never perfectly continuous.

Source: F. Galton (1886), "Regression towards mediocrity in hereditary stature", Journal of the Anthropological Institute 15, 246–263; data as transcribed by J. A. Hanley (2004), The American Statistician 58, 237–243.

From histogram to density curve

This is Example 5.1 in motion. We simulate adult men's heights (in inches) like the sons in Galton's data: mean 69.2, standard deviation 2.6. With few people and wide bins, the density histogram is blocky. Add more people and shrink the bins: the outline smooths into a curve. The probability of a bin is its area (height × width), so the probability of one exact height shrinks to zero as the bins shrink, even though the curve's height does not.

people = bin width = inches show the density curve

Probability density functions

The histogram shows the key idea. For a measurement, the probability of an interval is the area of the histogram over that interval. In the limit, the probability that $X$ falls between $a$ and $b$ is the area under the density curve from $a$ to $b$ .

Definition: Continuous Random Variable and Density

A random variable $X$ is continuous if there is a function $f (x)$ , called its probability density function (PDF, or just density), such that for all $a \leq b$ , $Pr (a \leq X \leq b) = \int_{a}^{b} f (x) d x .$ Every density satisfies

$f (x) \geq 0$ for all $x$ , and
$\int_{- \infty}^{\infty} f (x) d x$ $= 1$ : the total area is 1.

Three consequences follow straight away, and each is a common source of confusion.

Single values have probability zero. $Pr (X = a)$ $= \int_{a}^{a} f (x) d x$ $= 0$ . So for a continuous random variable, $Pr (a \leq X \leq b)$ $= Pr (a < X < b)$ : including the endpoints or not makes no difference, unlike the discrete case in Chapter 4.
A density is not a probability. $f (x)$ can be bigger than 1 (a density on an interval of length 0.1 must reach at least 10). What it measures is probability per unit of $x$ : for a tiny interval of width $Δ$ , $Pr (x \leq X \leq x$ $+ Δ)$ $\approx f (x) Δ$ . A taller curve means more probability packed into each inch.
Probability zero doesn't mean impossible. Every height has probability zero, yet one of them happens. (The dart in Section 3.3 made the same point.)

Example 5.2: Where Does the Dart Land? ★ lecture example

A dart lands at a uniformly random point on a round board of radius 1 (say, 1 foot). Let $R$ be its distance from the centre. (a) Show that $Pr (R \leq r)$ $= r^{2}$ for $0 \leq r \leq 1$ . (b) Find the density of $R$ and check that it's a valid density. (c) Find the probability that the dart lands within half a foot of the centre, and the probability it lands between 0.5 and 0.6 feet out.

(a) Probability is area. "Uniformly random point" means the probability of landing in any region is that region's share of the board's area. The dart lands within $r$ of the centre exactly when it lands in the circle of radius $r$ , so $Pr (R \leq r)$ $= \frac{area of circle of radius r}{area of the board}$ $= \frac{π r^{2}}{π \cdot 1^{2}}$ $= r^{2} .$

(b) The density is the rate at which this probability grows, its derivative: $f (r) = \frac{d}{d r} r^{2} = 2 r, 0 \leq r \leq 1,$ and $f (r) = 0$ outside $[0, 1]$ . It's never negative, and its total area is $\int_{0}^{1} 2 r d r$ $= [r^{2}]_{0}^{1} = 1$ . Valid.

(c) Within half a foot: $Pr (R \leq 0.5)$ $= {0.5}^{2}$ $= 0.25$ . Between 0.5 and 0.6: $Pr (0.5 \leq R \leq 0.6)$ $= \int_{0.5}^{0.6} 2 r d r$ $= {0.6}^{2}$ $- {0.5}^{2}$ $= 0.36$ $- 0.25$ $= 0.11 .$

The surprise. Only a quarter of random darts land within half the radius of the centre. The density $2 r$ rises towards the edge, because there is more board far out: a thin ring at distance $r$ has circumference $2 π r$ , so outer rings hold more area. Throw some darts:

set.seed(5)
x <- runif(100000, -1, 1); y <- runif(100000, -1, 1)
inside <- x^2 + y^2 <= 1                   # keep points that hit the board
R <- sqrt(x[inside]^2 + y[inside]^2)
mean(R <= 0.5)                             # close to 0.25
mean(R >= 0.5 & R <= 0.6)                  # close to 0.11

output

[1] 0.2482223
[1] 0.1093537

Picture. The density $f (r) = 2 r$ , with the two probabilities in (c) as areas:

r <- seq(0, 1, by = 0.005)
d <- data.frame(r, f = 2 * r)
ggplot(d, aes(r, f)) +
  geom_area(data = subset(d, r <= 0.5), fill = "steelblue", alpha = 0.4) +
  geom_area(data = subset(d, r >= 0.5 & r <= 0.6), fill = "orange", alpha = 0.7) +
  geom_line(linewidth = 1) +
  annotate("text", x = 0.36, y = 0.2, label = "0.25") +
  annotate("text", x = 0.55, y = 1.35, label = "0.11") +
  labs(x = "distance from centre, r (feet)", y = "density f(r)")

The CDF, again

The cumulative distribution function $F (x)$ $= Pr (X \leq x)$ is defined exactly as in Chapter 4. For a continuous random variable it is the area under the density to the left of $x$ : $F (x) = \int_{- \infty}^{x} f (t) d t .$ By the Fundamental Theorem of Calculus, the density is the derivative of the CDF, $f (x) = F' (x)$ , which is how we got $2 r$ from $r^{2}$ in Example 5.2. Instead of a staircase, the CDF of a continuous random variable is a smooth curve rising from 0 to 1, and it is steepest where the density is highest.

The CDF also gives quantiles. The $p$ -quantile of a continuous random variable is the value $q$ with $F (q) = p$ : a fraction $p$ of the probability lies to its left. The median is the 0.5-quantile.

Example 5.3: The Median Dart ★ lecture example

For the dart in Example 5.2, find the median distance from the centre and the 90th percentile.

Solve $F (q) = p$ with $F (r) = r^{2}$ : $q^{2} = p$ , so $q = \sqrt{p}$ .

Median: $q$ $= \sqrt{0.5}$ $\approx 0.707$ feet. Half of all darts land more than 0.7 feet from the centre of a 1-foot board.
90th percentile: $q$ $= \sqrt{0.9}$ $\approx 0.949$ feet. One dart in ten lands in the outermost 5% of the radius.

So a "typical" random dart is much nearer the rim than the centre, which is why darts players aim, and why a random scatter of shots on a target looks emptier in the middle than you'd expect.

Check Yourself

A continuous random variable has density

f

. Which statement is always true?

f (x) \leq 1

for every

x

Pr (X = 2)

= f (2)

Pr (1 \leq X \leq 2)

= Pr (1 < X < 2)

For what value of $c$ is $f (x) = c x^{2}$ on $0 \leq x \leq 1$ (and 0 elsewhere) a valid density?

c =

With that density, $f (x) = 3 x^{2}$ on $[0, 1]$ , find $Pr (X \leq \frac{1}{2})$ .

Pr(X ≤ 1/2) =

Think About It: How Precise Is a Measurement?

A forensic report says a footprint is "28.0 cm long". The probability that a true length is exactly 28.0 cm is zero, so what does this report actually tell a jury? How would you rewrite it to be honest about uncertainty?

It really means the length lies in an interval: perhaps 27.95 to 28.05 cm if it was rounded to the nearest millimetre, or wider still once the ruler's error and the stretch of the impression are counted. An honest report gives the interval or a measurement uncertainty: "28.0 cm ± 0.3 cm (95%)", say. This is exactly the shift this section makes: continuous quantities are described by probabilities of intervals, not of single values. Chapter 8 shows how to compute such intervals.

Key Takeaways

A continuous random variable has a density $f (x) \geq 0$ with total area 1; probabilities are areas: $Pr (a \leq X \leq b)$ $= \int_{a}^{b} f (x) d x$ .
Single values have probability zero, so $\leq$ and $<$ don't matter. A density is probability per unit, and can exceed 1.
$F (x)$ $= \int_{- \infty}^{x} f (t) d t$ and $f = F'$ . The $p$ -quantile solves $F (q) = p$ .

5.2 Expected value, variance and the uniform distribution

Every idea from Chapter 4 carries over, with sums replaced by integrals and the PMF replaced by the density.

Definition: Expected Value and Variance of a Continuous Random Variable

If $X$ has density $f$ , its expected value and variance are $E [X]$ $= \int_{- \infty}^{\infty} x f (x) d x,$ $Var [X]$ $= \int_{- \infty}^{\infty} {(x - μ)}^{2} f (x) d x,$ where $μ = E [X]$ (provided the integral of $| x | f (x)$ is finite). More generally, $E [g (X)]$ $= \int g (x) f (x) d x$ .

All the rules carry over unchanged, with the same proofs (sums become integrals): linearity $E [a X$ $+ b]$ $= a E [X]$ $+ b$ ; the shortcut $Var [X]$ $= E [X^{2}]$ $- (E [X])^{2}$ ; and $Var [a X$ $+ b]$ $= a^{2} Var [X]$ . The interpretations do too: $E [X]$ is the long-run average and the balance point of the density, and the standard deviation is a typical distance from the mean.

Example 5.4: The Average Dart ★ lecture example

For the dart distance $R$ with density $f (r) = 2 r$ on $[0, 1]$ , find $E [R]$ , $Var [R]$ and $SD [R]$ . Compare the mean with the median from Example 5.3.

Mean. $E [R]$ $= \int_{0}^{1} r$ $\cdot 2 r d r$ $= \int_{0}^{1} 2 r^{2} d r$ $= {[\frac{2 r^{3}}{3}]}_{0}^{1}$ $= \frac{2}{3}$ $\approx 0.667 .$

Variance, by the shortcut formula. First $E [R^{2}]$ $= \int_{0}^{1} r^{2}$ $\cdot 2 r d r$ $= \int_{0}^{1} 2 r^{3} d r$ $= {[\frac{r^{4}}{2}]}_{0}^{1}$ $= \frac{1}{2},$ so $Var [R]$ $= \frac{1}{2}$ $- {(\frac{2}{3})}^{2}$ $= \frac{1}{2}$ $- \frac{4}{9}$ $= \frac{1}{18},$ $SD [R]$ $= \sqrt{1 / 18}$ $\approx 0.236 .$

Mean vs median. The mean, 0.667, is below the median, 0.707. The density rises to the right and stops abruptly at 1, so its long tail is on the left: it's left-skewed, and (as for data in Chapter 2) the mean is pulled towards the long tail. R can check the integrals numerically:

integrate(function(r) r * 2 * r, lower = 0, upper = 1)      # E[R]
integrate(function(r) r^2 * 2 * r, lower = 0, upper = 1)    # E[R^2]

output

0.6666667 with absolute error < 7.4e-15
0.5 with absolute error < 5.6e-15

ggplot(data.frame(r = seq(0, 1, by = 0.005)), aes(r, 2 * r)) +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = 2/3) +
  geom_vline(xintercept = sqrt(0.5), linetype = "dashed") +
  annotate("text", x = 2/3 - 0.02, y = 1.9, label = "mean", hjust = 1) +
  annotate("text", x = sqrt(0.5) + 0.02, y = 1.9, label = "median", hjust = 0) +
  labs(x = "distance from centre, r (feet)", y = "density f(r)")

The uniform distribution

The simplest continuous distribution spreads probability evenly over an interval.

Definition: Uniform Distribution

$X$ has the uniform distribution on $[a, b]$ , written $X \sim Unif (a, b)$ , if its density is constant on that interval: $f (x) = \frac{1}{b - a}, a \leq x \leq b,$ and 0 elsewhere. Then $E [X] = \frac{a + b}{2}$ and $Var [X]$ $= \frac{{(b - a)}^{2}}{12}$ .

The total area is a rectangle of width $b - a$ and height $\frac{1}{b - a}$ , which is 1. For the mean, $\int_{a}^{b} \frac{x}{b - a} d x$ $= \frac{b^{2} - a^{2}}{2 (b - a)}$ $= \frac{a + b}{2}$ , the midpoint, as symmetry suggests. For the variance, $E [X^{2}]$ $= \int_{a}^{b} \frac{x^{2}}{b - a} d x$ $= \frac{b^{3} - a^{3}}{3 (b - a)}$ $= \frac{a^{2} + a b + b^{2}}{3}$ , and subtracting ${(\frac{a + b}{2})}^{2}$ leaves $\frac{{(b - a)}^{2}}{12}$ after a little algebra. $◼$ ∎

For a uniform distribution, probabilities are just lengths: $Pr (c \leq X \leq d)$ $= \frac{d - c}{b - a}$ for any $[c, d]$ inside $[a, b]$ .

Example 5.5: Waiting for the Bus ★ lecture example

A bus comes every 12 minutes, like clockwork, and you arrive at the stop at a random time without checking the timetable. Let $W$ be your wait. (a) What is the distribution of $W$ ? (b) What is your expected wait, and its standard deviation? (c) What is the probability you wait more than 10 minutes?

(a) Arriving "at a random time" makes every point in the 12-minute cycle equally likely, so $W$ $\sim Unif (0,$ $12)$ , with density $\frac{1}{12}$ on $[0, 12]$ .

(b) $E [W]$ $= \frac{0 + 12}{2}$ $= 6$ minutes, and $SD [W]$ $= \sqrt{\frac{12^{2}}{12}}$ $= \sqrt{12}$ $\approx 3.46$ minutes.

1 - punif(10, min = 0, max = 12)

output

[1] 0.1666667

Critical thinking. Real buses aren't like clockwork: they bunch up, and gaps vary. Surprisingly, irregular buses make your expected wait longer than half the average gap, because you're more likely to arrive during a long gap than a short one (long gaps take up more of the timeline). This is called the inspection paradox, and it's why waits always feel longer than the timetable suggests.

ggplot() +
  annotate("rect", xmin = 10, xmax = 12, ymin = 0, ymax = 1/12, fill = "orange", alpha = 0.7) +
  annotate("segment", x = 0, xend = 12, y = 1/12, yend = 1/12, linewidth = 1) +
  geom_vline(xintercept = 6, linetype = "dashed") +
  scale_y_continuous(limits = c(0, 0.12)) +
  labs(x = "wait W (minutes)", y = "density")

Simulation: turning uniform numbers into anything

Computers generate random numbers that are uniform on $[0, 1]$ ; runif() in R. How does R produce Gaussian or exponential random numbers from them? The trick uses quantiles: if $U \sim Unif (0, 1)$ and $F$ is a continuous CDF, then $X = F^{- 1} (U)$ , the $U$ -quantile, has CDF $F$ .

For any $x$ , $X \leq x$ exactly when $F^{- 1} (U) \leq x$ , that is, when $U \leq F (x)$ (since $F$ is increasing). And $Pr (U \leq F (x))$ $= F (x)$ because $U$ is uniform. $◼$ ∎

For the dart, $F (r) = r^{2}$ , so $R = \sqrt{U}$ simulates dart distances without throwing any darts. This is the inverse transform method.

Turning random numbers into any distribution

A computer's basic random numbers are uniform between 0 and 1. To get a draw from another distribution, pick a uniform number u on the vertical axis, go across to the CDF, and drop down to the horizontal axis: the x you land on is the u-quantile. Press draw to do this many times; the histogram of the results takes the shape of the target density, because steep parts of the CDF catch many u's.

target distribution

Think About It: The Inspection Paradox in the Wild

A university reports that its average class size is 30. A student survey finds that the average student is in a class of 60. Both are correct (Chapter 2's "pooling in your own life" box raised the same puzzle; now we can put numbers on it). How? Which number would you use to predict the size of a class a randomly chosen student walks into, and which to plan the number of classrooms?

Big classes contain many students, so when you sample students, you hit big classes more often, just as a random arrival at a bus stop is more likely to fall in a long gap. If there are nine classes of 20 and one of 120, the average class has $300 / 10 = 30$ students, but 120 of the 300 students are in the big class, so the average student's class has $(180 \times 20$ $+ 120 \times 120) / 300$ $= 60$ students. To predict a random student's experience, use the student-weighted average; to plan rooms, the class average. Both are predictions; they just predict different things. Before using any average, ask: averaged over what?

Key Takeaways

$E [X]$ $= \int x f (x) d x$ and $Var [X]$ $= E [X^{2}]$ $- μ^{2}$ : the same rules as Chapter 4, with integrals.
Uniform $Unif (a, b)$ : flat density; mean $\frac{a + b}{2}$ , variance $\frac{{(b - a)}^{2}}{12}$ ; probabilities are lengths.
If $U$ is uniform, $F^{- 1} (U)$ has CDF $F$ : that's how computers simulate any distribution.

Lecture 13 — The Gaussian distribution

5.3 The Gaussian distribution

The bell curve

Galton's sons' heights, measurement errors in Michelson's speed of light, birth weights, exam scores: an astonishing number of measurements have the same symmetric, bell-shaped distribution. It has a formula, and it is the most important distribution in statistics.

Definition: Gaussian Distribution

$X$ has the Gaussian distribution with mean $μ$ and standard deviation $σ > 0$ , written $X \sim N (μ, σ^{2})$ , if its density is $f (x)$ $= \frac{1}{σ \sqrt{2 π}} exp (- \frac{{(x - μ)}^{2}}{2 σ^{2}}),$ $- \infty$ $< x$ $< \infty .$ Then $E [X] = μ$ and $Var [X] = σ^{2}$ . The standard Gaussian distribution is $N (0, 1)$ , and a standard Gaussian random variable is usually called $Z$ .

(The notation follows the convention that the second number is the variance $σ^{2}$ . R's functions take the standard deviation instead: pnorm(x, mean, sd). Always check which one a formula or program wants.)

History: Why Gaussian and Not Normal?

Most textbooks, and R, call this distribution the normal distribution. We keep the standard notation $N (μ, σ^{2})$ , and R's functions dnorm, pnorm, qnorm and rnorm, where the "norm" stands for normal.

The word normal was applied to the curve by several scientists in the 1870s and 1880s, including Francis Galton, and Karl Pearson made it standard in the 1890s. It carried a message. In the 1830s and 1840s the Belgian astronomer Adolphe Quetelet had argued that the average of a human trait described an ideal "average man" (l'homme moyen), and that variation around it was like error in a measurement: nature missing its target. Galton used the bell curve to rank people into grades of "natural ability", and the language of a normal middle and abnormal tails fitted the eugenic project of sorting people into the fit and the unfit (the history box after Example 5.6 tells that story). Pearson himself came to regret the name. In 1920 he wrote:

"Many years ago I called the Laplace–Gaussian curve the normal curve, which name, while it avoids an international question of priority, has the disadvantage of leading people to believe that all other distributions of frequency are in one sense or another 'abnormal'."

The name Gaussian is already standard in physics, engineering and machine learning, and it describes the curve without judging anyone. It isn't a perfect name either. Abraham de Moivre found the curve in 1733, and Pierre-Simon Laplace developed its theory alongside Gauss. This is an example of what the statistician Stephen Stigler called his "law of eponymy": no scientific discovery is named after its original discoverer. But it is a neutral name, and it carries a useful reminder: most data are not Gaussian, and that isn't a defect.

When you read other books and papers, translate:

these notes	other books
Gaussian	normal
standard Gaussian	standard normal
Gaussian Q-Q plot	normal Q-Q plot
log-Gaussian	lognormal

How to read the formula: the curve is highest at $x = μ$ , where the exponent is zero. As $x$ moves away from $μ$ , the exponent becomes more and more negative quadratically, so the density falls off very fast: values more than 3 standard deviations from the mean are rare, and beyond 5 they are almost unheard of. The factor $\frac{1}{σ \sqrt{2 π}}$ in front makes the total area 1 (a famous integral, proved in the exercises of Chapter 6), and the $σ$ inside sets the width.

Standardizing

There is no formula for the area under a Gaussian curve; it has to be computed numerically. Fortunately, every Gaussian distribution is a stretched and shifted copy of the standard Gaussian, so one table (or one R function) covers them all.

Theorem: Standardizing a Gaussian Random Variable

If $X \sim N (μ, σ^{2})$ , then $Z$ $= \frac{X - μ}{σ}$ $\sim N (0,$ $1)$ . Consequently $Pr (X \leq x)$ $= Pr (Z \leq \frac{x - μ}{σ})$ $= Φ (\frac{x - μ}{σ}),$ where $Φ$ (capital phi) is the CDF of the standard Gaussian.

For any $z$ , $\begin{matrix} Pr (Z \leq z) & = Pr (X \leq μ + σ z) \\ = \int_{- \infty}^{μ + σ z} f (x) d x, \end{matrix}$ where $f$ is the $N (μ, σ^{2})$ density above. Substitute $t = (x - μ) / σ$ , so $d x = σ d t$ and the upper limit becomes $z$ : $Pr (Z \leq z) = \int_{- \infty}^{z} \frac{1}{\sqrt{2 π}} e^{- t^{2} / 2} d t,$ which is the standard Gaussian CDF. $◼$ ∎

The quantity $\frac{x - μ}{σ}$ is a z-score, just as for data in Chapter 2: the number of standard deviations $x$ lies from the mean. Standardizing turns every Gaussian question into a question about $Z$ .

Idea: The 68–95–99.7 Rule

For any Gaussian distribution, the probability of falling within 1, 2 and 3 standard deviations of the mean is about 68%, 95% and 99.7%: $Pr (| Z | \leq 1)$ $\approx 0.683,$ $Pr (| Z | \leq 2)$ $\approx 0.954,$ $Pr (| Z | \leq 3)$ $\approx 0.997 .$ This is where Chapter 2's rule of thumb comes from. The precise value for 95% is $\pm 1.96$ standard deviations, a number you'll meet again in Chapter 8.

z <- seq(-4, 4, by = 0.01); d <- data.frame(z, f = dnorm(z))
ggplot(d, aes(z, f)) +
  geom_area(data = subset(d, abs(z) <= 3), fill = "steelblue", alpha = 0.15) +
  geom_area(data = subset(d, abs(z) <= 2), fill = "steelblue", alpha = 0.25) +
  geom_area(data = subset(d, abs(z) <= 1), fill = "steelblue", alpha = 0.4) +
  geom_line(linewidth = 1) +
  annotate("text", x = 0, y = 0.17, label = "68.3%") +
  annotate("text", x = 2.45, y = 0.09, label = "95.4%", size = 3.8) +
  annotate("text", x = 3.35, y = 0.04, label = "99.7%", size = 3.8) +
  scale_x_continuous(breaks = -3:3, labels = c("μ−3σ", "μ−2σ", "μ−σ", "μ", "μ+σ", "μ+2σ", "μ+3σ")) +
  labs(x = NULL, y = "density")

In R, pnorm(x, mean, sd) gives $Pr (X \leq x)$ and qnorm(p, mean, sd) gives the $p$ -quantile.

Example 5.6: How Tall Were Galton's Sons? ★ lecture example

Galton's 481 sons had mean height 69.2 inches and standard deviation 2.6 inches. Model their heights as $N (69.2, {2.6}^{2})$ . (a) What fraction does the model predict between 67 and 72 inches? Compare with the data (Example 5.1). (b) What fraction were 6 ft 2 in (74 inches) or taller? (c) How tall must a son be to be in the tallest 5%? (d) Is the Gaussian model reasonable?

(a) Standardize both ends: $z_{1}$ $= \frac{67 - 69.2}{2.6}$ $\approx - 0.85,$ $z_{2}$ $= \frac{72 - 69.2}{2.6}$ $\approx 1.08,$ so $Pr (67 \leq X \leq 72)$ $= Φ (1.08)$ $- Φ (- 0.85)$ $\approx 0.860$ $- 0.198$ $= 0.662 .$ The data gave 0.657. Very close. (With the unrounded mean and SD, as in the R code below, the model gives 0.657 too.)

(b) $z$ $= \frac{74 - 69.2}{2.6}$ $\approx 1.85$ , so $Pr (X \geq 74)$ $= 1$ $- Φ (1.85)$ $\approx 0.032$ . In the data, 16 of 481 sons (3.3%) were at least 74 inches.

(c) We need the 95th percentile. The standard Gaussian's is $z = 1.645$ (qnorm(0.95)), so $x$ $= μ$ $+ z σ$ $= 69.2$ $+ 1.645$ $\times 2.6$ $\approx 73.5 inches .$

mu <- mean(sons); sigma <- sd(sons)
c(mu, sigma)
pnorm(72, mu, sigma) - pnorm(67, mu, sigma)      # (a) model
mean(sons >= 74); 1 - pnorm(74, mu, sigma)       # (b) data and model
qnorm(0.95, mu, sigma); quantile(sons, 0.95)     # (c) model and data

output

[1] 69.234096  2.623905
[1] 0.6568228
[1] 0.03326403
[1] 0.03465909
[1] 73.55004
95% 
 73

(d) The model's answers match the data closely, and a picture shows why:

ggplot(data.frame(sons), aes(sons)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 1, boundary = 0, colour = "white") +
  stat_function(fun = dnorm, args = list(mean = mu, sd = sigma), linewidth = 1) +
  labs(x = "height of son (inches)", y = "density")

A Gaussian Q-Q plot (Section 2.5), qqnorm(sons); qqline(sons), is nearly straight. The model describes the sons well. It would not describe all 934 children together, because sons and daughters are two groups with different means; mixing them gives a flatter, wider distribution than either.

Prediction and inference. The model was inferred from data (the mean and SD estimate the population's), and then used to predict fractions. Those predictions are checkable, and checking them against the data is how we decided the model was trustworthy. Its predictions far out in the tails, where there is little data to check against, deserve less trust.

History: Galton's Heights, and Why He Collected Them

The heights in this chapter come from a remarkable dataset. In 1884 Francis Galton (1822–1911), a half-cousin of Charles Darwin, collected the heights of the parents and adult children of 205 families, partly by offering cash prizes for detailed family records. From these data he found that tall parents tend to have tall children, but less extremely tall than themselves. He published this in 1886 as "Regression towards mediocrity in hereditary stature", and it is the origin of the words regression and regression to the mean (Section 4.7). Two years later he introduced correlation (Chapter 6). Every regression model fitted today, in medicine, economics or machine learning, descends from this work.

Galton collected the data for a reason. In 1883 he coined the word eugenics for a programme to improve the human population by encouraging people he considered "fit" to have more children and discouraging the "unfit". In Hereditary Genius (1869) he had ranked whole races by their supposed natural ability. Others turned the idea into policy:

more than 30 American states passed compulsory sterilization laws; the US Supreme Court upheld Virginia's law in Buck v. Bell (1927), and more than 60,000 people were sterilized under such laws;
eugenic arguments helped shape the US Immigration Act of 1924, which restricted immigration from southern and eastern Europe and barred it from Asia;
Nazi Germany's sterilization law of 1933 drew on American models.

In 2020 University College London, where Galton had endowed the chair of eugenics, removed his name and Pearson's from its buildings.

None of this followed from the mathematics. Regression and correlation measure how quantities vary together; they cannot say why. Galton read the similarity of parents and children as evidence that ability and character are mainly inherited. But families also share wealth, food, schooling and opportunity, which are confounders in the language of Chapter 1's DAGs, and his data could not separate these causes. The tools were sound; the inference was not.

The words carry the history too. "Regression" originally meant a drift back towards "mediocrity". That is one reason Aubrey Clayton, in Bernoulli's Fallacy (2021), proposes calling linear regression simply linear modelling, and avoiding "correlation" altogether. These notes keep the standard words, which you will need in order to read the scientific literature. But knowing their origin is part of being statistically literate. The same calculation can inform or mislead, and only a reader who understands where the numbers came from can tell which.

Causal Inference: Why Are So Many Things Gaussian?

Why should heights follow a bell curve? Think about what causes a person's height: hundreds of genes, each nudging height up or down a little, plus nutrition, childhood illnesses, sleep, and much else. Height is roughly the sum of many small, independent influences. The Central Limit Theorem of Chapter 7 says that such sums are approximately Gaussian, whatever the distributions of the individual pieces. That's the deep reason the bell curve is everywhere: measurement errors, test scores and birth weights are all sums of many small effects.

The same reasoning tells you when to expect a Gaussian distribution to fail. If one cause is large and dominant, the distribution follows that cause instead: adult heights for men and women together are a mixture of two bells, because sex has a big effect on height. Incomes, house prices and pollutant concentrations are shaped by causes that multiply rather than add, and are right-skewed (Section 5.6). Before assuming a Gaussian model, ask what produced the data: a DAG with many small arrows into the outcome suggests a bell; one big arrow suggests something else.

Example 5.7: Just Over the Limit? ★ lecture example

A breathalyser's readings follow a Gaussian distribution around the true blood alcohol concentration (BAC), with standard deviation 0.005 g/dL. A driver's true BAC is 0.075, below the legal limit of 0.08. (a) What is the probability that a single reading is 0.08 or above? (b) The policy is to take two independent readings and use the lower. What is the probability that both are 0.08 or above?

(a) The reading $X$ $\sim N (0.075,$ ${0.005}^{2})$ . Standardize: $Pr (X \geq 0.08)$ $= Pr (Z \geq \frac{0.08 - 0.075}{0.005})$ $= Pr (Z \geq 1)$ $= 1$ $- Φ (1)$ $\approx 0.159 .$ About one in six readings would put this driver over the limit, though their true level is under it.

(b) Both readings must be at least 0.08. They are independent, so multiply (Section 3.5): $Pr (both \geq 0.08)$ $= {0.159}^{2}$ $\approx 0.025 .$

1 - pnorm(0.08, mean = 0.075, sd = 0.005)        # (a)
(1 - pnorm(0.08, mean = 0.075, sd = 0.005))^2    # (b)

output

[1] 0.1586553
[1] 0.02517149

x <- seq(0.055, 0.095, by = 0.0001)
d <- data.frame(x, f = dnorm(x, 0.075, 0.005))
ggplot(d, aes(x, f)) +
  geom_area(data = subset(d, x >= 0.08), fill = "orange", alpha = 0.6) +
  geom_line() +
  geom_vline(xintercept = 0.08, linetype = "dashed") +
  labs(x = "breathalyser reading (g/dL)", y = "density")

What it means. Measurement error matters most for people near a threshold. A single reading would wrongly put a driver at 0.075 over the limit 16% of the time; requiring two readings cuts that to 2.5%. (The real situation is more complicated: the two readings aren't quite independent, because the same person and device produce both. Positive dependence would make the 2.5% an underestimate.) This is the kind of calculation forensic toxicologists use to set the "measurement uncertainty" they report alongside a result.

Turning it round. This example predicted readings from a known true level. A court faces the inference problem: given a reading of 0.081, what was the true level? Because readings scatter by 0.005, a true level anywhere from about 0.071 to 0.091 could easily produce it. Chapter 8 makes this precise with a confidence interval; a Bayesian would also use what is known about how common different BAC levels are among drivers who are tested.

The Gaussian approximation to the binomial

The binomial histograms in Chapter 4 looked more and more bell-shaped as $n$ grew. That's not a coincidence: a binomial count is a sum of $n$ independent indicators, so the Central Limit Theorem applies. For large $n$ , $Bin (n, p)$ is close to a Gaussian with the same mean $n p$ and variance $n p (1 - p)$ . Because the binomial is discrete, we treat each integer $k$ as the interval $[k - \frac{1}{2}, k + \frac{1}{2}]$ : this continuity correction noticeably improves the approximation.

Example 5.8: The Airline, One Last Time

The airline sells 110 tickets, each holder showing up with probability 0.9. Use the Gaussian approximation, with the continuity correction, to estimate the probability that more than 100 passengers show up, and compare with the exact binomial answer of 0.329 from Example 4.11.

$X$ $\sim Bin (110,$ $0.9)$ has $μ$ $= 110$ $\times 0.9$ $= 99$ and $σ$ $= \sqrt{110 \times 0.9 \times 0.1}$ $= \sqrt{9.9}$ $\approx 3.146$ . "More than 100" means $X \geq 101$ , which the continuity correction treats as "above 100.5": $Pr (X \geq 101)$ $\approx Pr (Z \geq \frac{100.5 - 99}{3.146})$ $= Pr (Z \geq 0.477)$ $\approx 0.317 .$

1 - pnorm(100.5, mean = 99, sd = sqrt(9.9))      # Gaussian approximation
1 - pbinom(100, 110, 0.9)                        # exact

output

[1] 0.3167767
[1] 0.3289856

The approximation gives 0.317 against the exact 0.329: close, though not perfect, because with $p = 0.9$ the binomial is slightly skewed. With a computer we'd just use the exact answer; the approximation earns its keep in Chapters 7 to 9, where it lets us reason about averages and proportions without knowing their exact distribution.

Picture. The Gaussian curve (with the same mean and SD) over the binomial bars; the bars for 101 or more passengers, and the area beyond 100.5, are what we compared:

k <- 88:110
d <- data.frame(k, p = dbinom(k, 110, 0.9), over = k >= 101)
ggplot(d, aes(k, p, fill = over)) + geom_col(width = 0.8) +
  stat_function(fun = dnorm, args = list(mean = 99, sd = sqrt(9.9)), inherit.aes = FALSE, linewidth = 1) +
  scale_fill_manual(values = c("grey60", "orange"), guide = "none") +
  labs(x = "passengers who show up", y = "probability")

History: Ceres, Gauss and the Birth of the Bell Curve

On 1 January 1801, the Italian astronomer Giuseppe Piazzi discovered a new object in the night sky: Ceres, now known as a dwarf planet. He tracked it for about six weeks, and then it disappeared into the glare of the Sun. Astronomers had only a few, slightly inaccurate measurements, and no idea where it would reappear.

A 24-year-old German mathematician, Carl Friedrich Gauss, took on the problem. Using new methods, including what became the method of least squares, he found the orbit that best fitted the inaccurate observations. In December 1801 astronomers pointed their telescopes where Gauss said Ceres would be, and there it was.

In 1809 Gauss explained why least squares works. He assumed that the average of repeated measurements is the most probable value of the true quantity, showed that this forces measurement errors to follow a particular symmetric, bell-shaped curve, and showed that least squares then gives the most probable orbit (Chapter 8's maximum likelihood). That is why the curve carries his name, though Abraham de Moivre had already found it in 1733, as an approximation to the binomial.

The bell curve began as a model of errors in measuring the sky. Applying it to people came later, with Quetelet, who measured the chests of Scottish soldiers, and then Galton. Both moves are described in the history boxes earlier in this section.

Check Yourself

Find $Pr (Z > 2)$ for a standard Gaussian $Z$ . (Use R: 1 - pnorm(2).)

Pr(Z > 2) =

Galton's daughters' heights are roughly $N (64.1, {2.4}^{2})$ inches. How tall would a daughter need to be to be taller than 99.99% of daughters? Give your answer to one decimal place.

height =

Test scores are

N (500, 100^{2})

. Which is more unusual: a score of 720, or a score of 290? 720290they are equally unusual

Think About It: Abnormal, or Just Unlucky?

A lab uses a "reference range" for a blood marker: the middle 95% of healthy adults, mean ± 2 SD. If a healthy patient has 20 independent tests, each with its own 95% range, what's the chance at least one result is flagged as abnormal? What does this mean for interpreting a single flagged result?

Each healthy result is outside its range with probability 0.05, so $Pr (at least one flag)$ $= 1$ $- {0.95}^{20}$ $\approx 0.64$ . Most perfectly healthy people having a 20-test panel will have something flagged. A single flag is weak evidence on its own: it's the prediction "95% of healthy people fall here" applied twenty times, so some misses are expected. This is the multiple-testing problem again (Paul the octopus, cancer clusters), and Chapter 9 returns to it. Doctors interpret a flag alongside symptoms and repeat tests, which is Bayes' theorem in practice.

Key Takeaways

Gaussian $N (μ, σ^{2})$ : symmetric bell, mean $μ$ , SD $σ$ . About 68%, 95% and 99.7% lies within 1, 2 and 3 SD.
Standardize: $Z$ $= (X$ $- μ) / σ$ $\sim N (0,$ $1)$ . In R: pnorm for probabilities, qnorm for quantiles.
Sums of many small independent effects are approximately Gaussian (the Central Limit Theorem). One dominant cause, or multiplying causes, break the bell.
$Bin (n,$ $p)$ $\approx N (n p,$ $n p (1$ $- p))$ for large $n$ ; use the continuity correction.

Lecture 14 — Exponential, gamma, beta and other shapes

5.4 Waiting times: exponential and gamma

The exponential distribution

In Chapter 4, the Poisson distribution counted rare events (calls, accidents, decays) in a window of time. Now ask a different question about the same events: how long until the next one? The answer is a continuous random variable, and if events occur independently at a constant rate $λ$ , its distribution is exponential.

Definition: Exponential Distribution

$X$ has the exponential distribution with rate $λ > 0$ , written $X \sim Exp (λ)$ , if its density is $f (x) = λ e^{- λ x}, x \geq 0,$ and 0 for $x < 0$ . Its CDF is $F (x)$ $= 1$ $- e^{- λ x}$ , so $Pr (X > x)$ $= e^{- λ x}$ , and $E [X]$ $= \frac{1}{λ},$ $Var [X]$ $= \frac{1}{λ^{2}},$ $median$ $= \frac{ln 2}{λ} .$

You've already met this density: Example 0.4 showed its total area is 1 and Example 0.5 computed its mean, $1 / λ$ . The CDF comes from integrating: $\int_{0}^{x} λ e^{- λ t} d t$ $= 1$ $- e^{- λ x}$ . Setting $F (m) = \frac{1}{2}$ gives the median $m$ $= ln 2 / λ$ $\approx 0.69 / λ$ , below the mean: the exponential is right-skewed, with many short waits and a few very long ones.

Theorem: Waiting Times in a Poisson Process

If events occur as a Poisson process with rate $λ$ (Section 4.7), the waiting time $T$ until the first event is $Exp (λ)$ .

The wait is longer than $t$ exactly when there are no events in $[0, t]$ . The number of events in that window is $Pois (λ t)$ , so $Pr (T > t)$ $= Pr (no events in [0, t])$ $= e^{- λ t} \frac{{(λ t)}^{0}}{0!}$ $= e^{- λ t},$ which is the exponential tail probability. $◼$ ∎

So the Poisson and exponential distributions are two views of the same process: the Poisson counts events in a fixed time, the exponential measures time between events. With 911 calls at 4.5 per hour (Example 4.20), the gap between calls is $Exp (4.5)$ , with mean $\frac{1}{4.5}$ hour $\approx 13.3$ minutes, and $Pr (gap > 20 min)$ $= e^{- 4.5 / 3}$ $= e^{- 1.5}$ $\approx 0.223$ , exactly the probability of no calls in 20 minutes that we found there.

Theorem: The Exponential Is Memoryless

If $X \sim Exp (λ)$ , then for all $s, t \geq 0$ , $Pr (X > s$ $+ t ∣ X > s)$ $= Pr (X > t) .$

${X > s + t}$ is contained in ${X > s}$ , so $Pr (X > s$ $+ t ∣ X > s)$ $= \frac{Pr (X > s + t)}{Pr (X > s)}$ $= \frac{e^{- λ (s + t)}}{e^{- λ s}}$ $= e^{- λ t}$ $= Pr (X > t) .$ $◼$ ∎

As for the geometric distribution (Section 4.5), past waiting tells you nothing about future waiting. If you've waited 30 minutes for a call, the remaining wait has the same distribution as if you'd just started. That's a strong assumption: it fits radioactive atoms and random arrivals, but not light bulbs or people, which wear out.

Example 5.9: British Coal-Mining Disasters ★ lecture example

R's boot package (installed with R) contains the dates of 191 explosions in British coal mines that killed 10 or more people, from 1851 to 1962. (a) For the disasters up to 1890, compute the gaps between consecutive disasters, in days, and check whether an exponential model fits. (b) Using the model, find the probability of a whole year without a disaster, and compare with the data. (c) What happened to the gaps after 1900?

(a) coal$date records each disaster as a decimal year. diff() gives the gaps between consecutive dates; multiplying by 365.25 converts to days.

library(boot)
gaps <- diff(coal$date) * 365.25            # days between disasters
early <- gaps[coal$date[-1] < 1890]         # gaps ending before 1890
length(early)
c(mean = mean(early), sd = sd(early), median = median(early))

output

[1] 122
    mean       sd   median 
115.5328 122.9763  80.5000

Three quick checks for an exponential distribution, all passed:

The mean (115.5 days) and standard deviation (123 days) should be equal, since both are $1 / λ$ . They nearly are.
The median should be $ln 2 \approx 0.69$ times the mean: $0.69$ $\times 115.5$ $\approx 80$ days. The data's median is 80.5.
The histogram should follow the curve $λ e^{- λ x}$ with $λ = 1 / 115.5$ :

ggplot(data.frame(days = early), aes(days)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 50, boundary = 0, colour = "white") +
  stat_function(fun = dexp, args = list(rate = 1 / mean(early)), linewidth = 1) +
  labs(x = "days between disasters (1851 to 1890)", y = "density")

(b) With $λ = 1 / 115.5$ per day, $Pr (gap > 365)$ $= e^{- 365 / 115.5}$ $\approx 0.042$ . In the data, 5 of the 122 gaps (4.1%) were longer than a year.

exp(-365 / mean(early))     # model
mean(early > 365)           # data
mean(gaps[coal$date[-1] >= 1900])   # (c) mean gap after 1900

output

[1] 0.04245645
[1] 0.04098361
[1] 408.2321

(c) After 1900, the average gap between disasters was about 408 days, more than three times as long. The rate of disasters fell sharply. A constant-rate exponential model fits each period, but not the whole century: the mechanism changed. Historians link the change to safety legislation, better ventilation and the decline of the most dangerous practices. (Whether the laws caused the drop is a causal question that needs more than these data: Chapter 1's DAGs again.)

Data: R. G. Jarrett (1979), "A note on the intervals between coal-mining disasters", Biometrika 66, 191–193.

Prediction and inference. Parts (a) and (b) are inference about the mechanism (disasters at a steady random rate, about one every 115 days) and a prediction from it (a disaster-free year about 4% of the time). Part (c) is the warning: a model fitted to 1851–1890 would have predicted the 1900s badly. Every prediction assumes the mechanism stays the same.

Example 5.10: How Old Is the Iceman? ★ lecture example

In 1991, hikers in the Alps found a frozen human body, now known as Ötzi the Iceman. Living things absorb carbon-14, a radioactive form of carbon, until they die; afterwards it decays. Each carbon-14 atom's lifetime is exponential, with a half-life of 5,730 years: the median lifetime. (a) Find the decay rate $λ$ . (b) What fraction of the original carbon-14 remains after $t$ years? (c) Ötzi's tissue contained about 53% of the carbon-14 found in living tissue. Estimate how long ago he died.

(a) The half-life is the median, $\frac{ln 2}{λ} = 5730$ , so $λ = \frac{ln 2}{5730} \approx 0.000121 per year .$ (The mean lifetime of an atom is $1 / λ \approx 8,267$ years, longer than the half-life, as for any exponential.)

(b) Each atom survives $t$ years with probability $Pr (X > t)$ $= e^{- λ t}$ . With trillions of atoms, the Law of Large Numbers makes the fraction surviving essentially equal to this probability: $fraction remaining$ $= e^{- λ t}$ $= {(\frac{1}{2})}^{t / 5730} .$ After one half-life, half remains; after two, a quarter; after three, an eighth.

(c) Solve $e^{- λ t} = 0.53$ : take logs, $- λ t = ln 0.53$ , so $t$ $= \frac{- ln 0.53}{λ}$ $= \frac{0.635}{0.000121}$ $\approx 5,250 years .$

lambda <- log(2) / 5730
lambda
-log(0.53) / lambda          # years since death

output

[1] 0.0001209681
[1] 5248.312

Ötzi died around 5,200 to 5,300 years ago, in the Copper Age, consistent with the radiocarbon dates that archaeologists published. (He had an arrowhead lodged in his shoulder: the oldest murder victim to be studied with modern forensic science.) The widget below shows the same process with 400 atoms: individual decays are completely unpredictable, but the fraction remaining is not.

t <- seq(0, 15000, by = 50)
ggplot(data.frame(t, frac = exp(-lambda * t)), aes(t, frac)) +
  geom_line(linewidth = 1) +
  annotate("segment", x = 0, xend = -log(0.53) / lambda, y = 0.53, yend = 0.53, linetype = "dashed") +
  annotate("segment", x = -log(0.53) / lambda, xend = -log(0.53) / lambda, y = 0, yend = 0.53, linetype = "dashed") +
  labs(x = "years since death", y = "fraction remaining")

This is inference in its purest form: an unobservable quantity (the time since death) inferred from a measurement (the fraction remaining) through a probability model of the mechanism (exponential decay). Real radiocarbon dates also report an uncertainty, because the 53% is itself measured with error.

Radioactive decay, one atom at a time

Each of the 400 squares is a radioactive atom. Every atom's lifetime is exponential with the same half-life, independently of the others: no atom "knows" how old it is. Press let time pass and watch the atoms decay (they turn grey). The right panel tracks the fraction still undecayed against the exponential curve e^−λt. Notice that the fraction halves in every half-life, whether it started at 100% or 25%.

The gamma distribution

The exponential is the wait for the first event. The wait for the $r$ th event is a sum of $r$ independent exponential gaps, and has a gamma distribution, just as the negative binomial was a sum of geometric waits.

Definition: Gamma Distribution

$X$ has the gamma distribution with shape $α > 0$ and rate $λ > 0$ , written $X$ $\sim Gamma (α,$ $λ)$ , if $f (x)$ $= \frac{λ^{α}}{Γ (α)} x^{α - 1} e^{- λ x},$ $x$ $> 0,$ where $Γ (α)$ $= \int_{0}^{\infty} t^{α - 1} e^{- t} d t$ is the gamma function ( $Γ (r) = (r - 1)!$ for whole numbers $r$ ). Then $E [X] = α / λ$ and $Var [X]$ $= α / λ^{2}$ . The exponential is the case $α = 1$ ; for a whole number $α = r$ , $X$ is the waiting time for the $r$ th event of a Poisson process with rate $λ$ .

The gamma is also a flexible model for any positive, right-skewed measurement (rainfall totals, insurance claims, hospital stays), whether or not it's literally a waiting time.

Example 5.11: The Third Call

911 calls arrive at 4.5 per hour. What is the probability that the third call of the afternoon arrives within 20 minutes? Answer it two ways.

As a waiting time. The wait $T$ for the 3rd call is $Gamma (3,$ $4.5)$ (in hours), and we want $Pr (T \leq \frac{1}{3})$ : pgamma(1/3, shape = 3, rate = 4.5).

As a count. The 3rd call arrives within 20 minutes exactly when there are at least 3 calls in those 20 minutes. The count is $Pois (4.5 / 3)$ $= Pois (1.5)$ , so $Pr (T \leq \frac{1}{3})$ $= 1$ $- Pr (N \leq 2)$ $= 1$ $- e^{- 1.5} (1 + 1.5 + \frac{{1.5}^{2}}{2})$ $\approx 1$ $- 0.809$ $= 0.191 .$

pgamma(1/3, shape = 3, rate = 4.5)
1 - ppois(2, 1.5)

output

[1] 0.1911532
[1] 0.1911532

Both give 0.191. Counting events and timing them are two descriptions of one process; whichever is easier to compute will do.

Check Yourself

Light bulbs in a hallway fail at random, with lifetimes exponential with mean 10 months. What is the probability a bulb lasts more than 10 months?

Pr(X > 10) =

The same bulb has already lasted 10 months. What is the probability it lasts at least another 10?

Pr =

Think About It: Memoryless or Not?

Which of these are plausibly memoryless (exponential), and which are not? (a) The time until a uranium atom decays. (b) A person's remaining lifetime. (c) The time until a car breaks down. (d) The time until the next customer enters a busy shop. For those that aren't, does the chance of the event in the next month go up or down with age?

(a) Memoryless: atoms have no internal clock, which is exactly what makes radiocarbon dating work. (b) Not: the chance of dying in the next year rises with age (after childhood), so lifetimes are better described by a Weibull or similar model with an increasing failure rate. (c) Mostly not: parts wear, so risk rises with age, though a new car may also have early "teething" failures. (d) Roughly memoryless over short periods, but the rate changes through the day, so a model with a constant $λ$ only works for a short window. Choosing a model is a claim about the mechanism; the exponential claims there is no ageing at all.

Key Takeaways

Exponential $Exp (λ)$ : the wait for the first event of a Poisson process. $Pr (X > x)$ $= e^{- λ x}$ ; mean $1 / λ$ ; median $ln 2 / λ$ .
Memoryless: no wear and tear. Half-lives are exponential medians.
Gamma $Gamma (α, λ)$ : the wait for the $α$ th event; mean $α / λ$ . Counting events and timing them give the same answers.

5.5 The beta distribution: uncertain proportions

Some quantities are proportions: a player's true free-throw percentage, the fraction of a drug seizure that's genuine, the share of voters supporting a policy. They live between 0 and 1, and the beta distribution is the standard model for them.

Definition: Beta Distribution

$X$ has the beta distribution with shape parameters $α, β > 0$ , written $X$ $\sim Beta (α,$ $β)$ , if $f (x)$ $= \frac{Γ (α + β)}{Γ (α) Γ (β)} x^{α - 1} {(1 - x)}^{β - 1},$ $0$ $< x$ $< 1 .$ Its mean is $E [X] = \frac{α}{α + β}$ and its variance $\frac{α β}{{(α + β)}^{2} (α + β + 1)}$ . $Beta (1, 1)$ is the uniform distribution on $[0, 1]$ .

The shape is very flexible: symmetric when $α = β$ , piled up near 1 when $α > β$ , near 0 when $α < β$ , and more concentrated as $α + β$ grows. Explore it:

The beta distribution's most important job is in Bayesian statistics, as a description of what we know about an unknown proportion.

Bayesian Inference: Learning a Free-Throw Percentage

In Chapter 4's Bayesian box, we learned a success probability $p$ from a grid of 11 candidate values. The beta distribution lets $p$ be any value between 0 and 1. And it has a remarkable property: if your prior for $p$ is $Beta (α, β)$ and you observe $s$ successes and $f$ failures in independent trials, your posterior is $Beta (α + s, β + f) .$ The prior's parameters act like "imaginary" earlier successes and failures, and data simply add to them. (The proof multiplies the prior density by the binomial likelihood $p^{s} {(1 - p)}^{f}$ : the result has the beta form, with the exponents increased by $s$ and $f$ .)

Example. A new basketball player makes 7 of her first 10 free throws. Starting from a uniform prior, $Beta (1, 1)$ , her posterior is $Beta (8, 4)$ :

p <- seq(0, 1, by = 0.005)
d <- rbind(data.frame(p, density = dbeta(p, 1, 1), which = "prior Beta(1, 1)"),
           data.frame(p, density = dbeta(p, 8, 4), which = "posterior Beta(8, 4)"))
ggplot(d, aes(p, density, linetype = which)) + geom_line(linewidth = 1) +
  labs(x = "free-throw success probability p", y = "density", linetype = NULL)
8 / (8 + 4)                         # posterior mean
1 - pbeta(0.5, 8, 4)                # Pr(p > 0.5 | data)
qbeta(c(0.025, 0.975), 8, 4)        # middle 95% of the posterior

output

[1] 0.6666667
[1] 0.8867188
[1] 0.3902574 0.8907366

Her posterior mean is $\frac{8}{12} \approx 0.67$ , a little below her observed 70%, pulled slightly towards the prior's 50%. There's an 89% probability that she's better than a 50% shooter, and her true percentage is 95% likely to lie between 0.39 and 0.89: after only 10 shots, we know surprisingly little. After 100 shots the posterior would be much narrower. This "credible interval" is the Bayesian cousin of Chapter 8's confidence interval.

Think About It: What Does the Prior Do?

Two scouts watch the same new player make 7 of 10 free throws. One starts from a uniform prior, $Beta (1, 1)$ ; the other, who knows that most professional players make 70% to 80%, starts from $Beta (15, 5)$ . Find each posterior mean. Whose estimate would you trust more after 10 shots, and after 1,000?

Posteriors: $Beta (8, 4)$ with mean $8 / 12 \approx 0.67$ , and $Beta (22, 8)$ with mean $22 / 30 \approx 0.73$ . After 10 shots the second scout's estimate is probably better, because her prior contains real information (the equivalent of 20 earlier shots) and 10 shots are very few. After 1,000 shots, the data overwhelm both priors: with, say, 700 makes the posteriors are $Beta (701, 301)$ and $Beta (715, 305)$ , with means 0.700 and 0.701. A prior matters when data are scarce, and it should be based on evidence that can be stated and checked; with plenty of data, reasonable priors agree.

5.6 More shapes of chance

Log-Gaussian: when causes multiply

Some quantities grow by percentages: a river gathers tributaries that each enlarge it, a city grows by a few percent a year, a pollutant's concentration is diluted by successive factors. When many independent effects multiply, the logarithm of the result is a sum of many effects, and so (by the reasoning of the causal box in Section 5.3) approximately Gaussian.

Definition: log-Gaussian Distribution

$X$ has the log-Gaussian distribution if $ln X$ $\sim N (μ,$ $σ^{2})$ . Then $X$ is positive and right-skewed, with median $e^{μ}$ and mean $e^{μ + σ^{2} / 2}$ , which is larger than the median. In R: plnorm(x, meanlog = mu, sdlog = sigma).

Example 5.12: How Long Is a River? ★ lecture example

R's rivers dataset gives the lengths, in miles, of 141 major rivers in North America (from the World Almanac, 1975). (a) Describe the distribution of lengths, and of their logarithms. (b) Fit a log-Gaussian model and use it to estimate the fraction of rivers longer than 1,000 miles. Compare with the data.

(a) The lengths are strongly right-skewed: most are a few hundred miles, but the Mississippi, Missouri and a few others run to thousands. On a log scale the skew largely disappears:

d <- rbind(data.frame(scale = "length (miles)", value = rivers),
           data.frame(scale = "log of length", value = log(rivers)))
ggplot(d, aes(value)) + geom_histogram(bins = 20, colour = "white") +
  facet_wrap(~ scale, scales = "free") + labs(x = NULL, y = "number of rivers")

(b) Fit a Gaussian to the logs: their mean is 6.18 and SD 0.59. Then $Pr (X > 1000)$ $= Pr (ln X > ln 1000)$ $= 1$ $- Φ (\frac{6.908 - 6.176}{0.591})$ $= 1$ $- Φ (1.24)$ $\approx 0.108 .$

m <- mean(log(rivers)); s <- sd(log(rivers))
c(m, s)
1 - plnorm(1000, meanlog = m, sdlog = s)   # model
mean(rivers > 1000)                         # data: 16 of 141
exp(m); median(rivers)                      # model median vs data median

output

[1] 6.1758789 0.5914841
[1] 0.1079775
[1] 0.1134752
[1] 481.0056
[1] 425

The model gives 10.8%; the data, 11.3% (16 rivers). The model's median, $e^{6.176} \approx 481$ miles, is a bit above the data's 425, a sign that the fit isn't perfect, but it captures the long tail far better than a Gaussian model could: a Gaussian with the data's mean (591) and SD (494) would give a 12% chance of a river with negative length.

In toxicology, drug and pollutant concentrations in the body are usually log-Gaussian, which is why toxicologists report geometric means and work on a log scale; ozone levels (Chapter 2) behave the same way.

Weibull: things that wear out

The exponential distribution says that a component's chance of failing next month doesn't depend on its age. Most things aren't like that: tyres, bones and hard drives wear out, so their failure rate rises with age; and some products fail early from manufacturing defects, so their failure rate falls. The Weibull distribution, with $Pr (X > x)$ $= e^{- {(x / θ)}^{k}}$ , handles both: shape $k > 1$ means wearing out, $k < 1$ means early failures, and $k = 1$ is the exponential. Engineers use it for reliability, and medical researchers for survival times.

Chi-squared and t: distributions for later

Two more distributions arise not as models of the world but from calculations on Gaussian data, and they'll be essential in Chapters 8 and 9.

Chi-squared with $k$ degrees of freedom: the distribution of $Z_{1}^{2} + \dots + Z_{k}^{2}$ , a sum of squares of $k$ independent standard Gaussians. It describes sample variances. (It's a gamma distribution with $α = k / 2$ and $λ = \frac{1}{2}$ .)
Student's t with $k$ degrees of freedom: like a standard Gaussian, but with heavier tails, reflecting extra uncertainty when the standard deviation must be estimated from a small sample. As $k$ grows it approaches the Gaussian.

History: Student, Who Worked at a Brewery

The t distribution was discovered by William Sealy Gosset, a chemist at the Guinness brewery in Dublin. Guinness needed to judge the quality of barley and hops from very small samples, and Gosset worked out how much the Gaussian distribution's predictions had to be widened when the standard deviation was estimated from just a handful of measurements. Guinness didn't allow its scientists to publish under their own names (after an earlier employee's paper had given away trade secrets, the company banned such publications), so his 1908 paper appeared under the pen name "Student", and the distribution has been "Student's t" ever since.

Think About It: Which Scale Should We Predict on?

A toxicologist reports that the average concentration of a pollutant in 50 water samples is 12 µg/L, but the median is 5 µg/L and the legal limit is 10. A water company says "typical water is well under the limit"; an activist says "on average the water is over the limit". Who is right? What summary would you report, and what would you want to predict?

Both are describing a right-skewed (probably log-Gaussian) distribution truthfully: most samples are low, but a few very high ones pull the mean up. Which summary matters depends on the question. For a person's long-run exposure, the mean is relevant, since high and low days average out over time. For the risk on a given day, we'd want the probability that a sample exceeds the limit, which a log-Gaussian model fitted to the logged data can estimate (as in Example 5.12). Report the median, the mean and the fraction over the limit, and show a histogram on a log scale. Neither single number tells the whole story.

5.7 Choosing a continuous model

Each named distribution has a characteristic shape. Seeing them side by side is the quickest way to match a model to a histogram. (The examples shown are a gamma with shape 3, a beta with parameters 2 and 5, a log-Gaussian with $σ = 0.6$ and a Weibull with shape 2; changing the parameters changes the details, not the family resemblance.)

curves <- function(name, x, y) data.frame(model = name, x = x, y = y)
d <- rbind(
  curves("uniform", seq(-0.2, 1.2, by = 0.005), dunif(seq(-0.2, 1.2, by = 0.005))),
  curves("Gaussian", seq(-4, 4, by = 0.02), dnorm(seq(-4, 4, by = 0.02))),
  curves("exponential", seq(0, 5, by = 0.02), dexp(seq(0, 5, by = 0.02))),
  curves("gamma", seq(0, 10, by = 0.02), dgamma(seq(0, 10, by = 0.02), 3)),
  curves("beta", seq(0, 1, by = 0.005), dbeta(seq(0, 1, by = 0.005), 2, 5)),
  curves("log-Gaussian", seq(0, 5, by = 0.02), dlnorm(seq(0, 5, by = 0.02), 0, 0.6)),
  curves("Weibull", seq(0, 3, by = 0.01), dweibull(seq(0, 3, by = 0.01), 2)))
d$model <- factor(d$model, levels = unique(d$model))
ggplot(d, aes(x, y)) + geom_line(linewidth = 1) +
  facet_wrap(~ model, scales = "free", ncol = 3) +
  scale_x_continuous(n.breaks = 3) + scale_y_continuous(n.breaks = 3) +
  labs(x = NULL, y = "density")

distribution	the story	values	mean
$Unif (a, b)$	equally likely anywhere in an interval	$[a, b]$	$\frac{a + b}{2}$
$N (μ, σ^{2})$	sums of many small independent effects; measurement error	all real numbers	$μ$
$Exp (λ)$	wait for the first random event; no ageing	$x \geq 0$	$\frac{1}{λ}$
$Gamma (α, λ)$	wait for the $α$ th event; positive skewed amounts	$x > 0$	$\frac{α}{λ}$
$Beta (α, β)$	an uncertain proportion	$(0, 1)$	$\frac{α}{α + β}$
log-Gaussian	products of many effects; concentrations, sizes	$x > 0$	$e^{μ + σ^{2} / 2}$
Weibull	lifetimes with wear-out or early failure	$x \geq 0$

Every distribution has d, p, q and r functions in R, exactly as in Chapter 4, named after its short name: unif, norm, exp, gamma, beta, lnorm and weibull (for example dexp, pexp, qexp, rexp).

Idea: Algorithm: Which Continuous Model?

What values are possible? All real numbers (Gaussian, t); positive only (exponential, gamma, log-Gaussian, Weibull); between 0 and 1 (beta); a fixed interval (uniform).
What mechanism produced it? Many small additive effects: Gaussian. Multiplicative effects: log-Gaussian. Waiting for random events: exponential or gamma. Wear-out: Weibull.
Look at the data: a histogram, and a Q-Q plot against the candidate model.
Check the claims you care about. A model can fit the middle and fail in the tails, which is often exactly where the important questions are.

Check Yourself

The time between successive cyber-attacks on a company's network, if attacks arrive independently at a steady average rate.

The unknown proportion of a large drug seizure that is genuine, as described by a Bayesian analyst after testing some bags.

The concentration of a pollutant in water samples, which is positive and strongly right-skewed, with a roughly bell-shaped histogram of its logarithm.

Prediction and inference in this chapter

Prediction. A fitted density predicts how often future values fall in any range: the fraction of sons over 74 inches, the chance of a year without a mine disaster, the probability that an honest driver reads over the limit, the share of rivers longer than 1,000 miles. Each time, we checked the prediction against data, and the agreement was the reason to trust the model.
Inference. Measurements are clues to hidden quantities. A carbon-14 fraction gave Ötzi's age; a breathalyser reading is evidence about a true alcohol level; seven free throws in ten updated our beliefs about a shooting percentage (the Bayesian route, with the beta distribution). The coal-mine gaps showed inference about change: the rate of disasters fell after 1900, and a model fitted to the early years would have predicted badly later.
The common thread. In every case the model's parameters (a mean, an SD, a rate, a proportion) were estimated from data, and we treated the estimates as if they were exact. They aren't: another 481 sons would give a slightly different mean. Chapter 6 develops the tools for several random variables at once, and Chapter 7 uses them to measure how much an estimate varies from sample to sample, which is the key to Chapters 8 and 9.

5.8 Interludes

Interlude (Optional): The Cauchy Distribution: When Averages Fail

A lighthouse stands 1 mile off a straight coast and rotates at a steady speed, sending out brief flashes at uniformly random angles. Where along the coast do the flashes hit? If the angle is $Θ$ $\sim Unif (- \frac{π}{2}, \frac{π}{2})$ , the flash hits the coast at $X = tan Θ$ . This $X$ has the Cauchy distribution, with density $f (x)$ $= \frac{1}{π (1 + x^{2})},$ $- \infty$ $< x$ $< \infty .$ It looks like a bell curve centred at 0, a little lower and wider. But its tails fall off so slowly (like $1 / x^{2}$ , not like $e^{- x^{2}}$ ) that the integral $\int | x | f (x) d x$ is infinite: the Cauchy distribution has no expected value. When the beam is nearly parallel to the coast, the flash lands enormously far away, and those rare enormous values dominate any average.

The consequence is startling: the Law of Large Numbers fails. The average of $n$ Cauchy draws has exactly the same Cauchy distribution as a single draw, however large $n$ is. Averaging doesn't help at all.

Why it matters: the Law of Large Numbers and the Central Limit Theorem need a finite mean (and variance). Real data with extremely heavy tails (insurance losses from catastrophes, the sizes of internet traffic bursts, losses in financial crashes) can behave more like the Cauchy than the Gaussian, and averages of such data are much less trustworthy than the Gaussian theory suggests.

Interlude (Optional): Families of Distributions

Many of this chapter's distributions come in families built from one basic shape.

Location–scale families. If $Z$ has some standard distribution, then $X = μ + σ Z$ has the same shape, shifted by $μ$ and stretched by $σ$ . The Gaussian family is one ( $Z \sim N (0, 1)$ ); so are the uniform, Cauchy and t families. Standardizing ( $Z$ $= (X$ $- μ) / σ$ ) works for any location–scale family, and probabilities for every member come from one table for $Z$ .
Exponential families. Most of the named distributions in Chapters 4 and 5 (binomial, Poisson, geometric, Gaussian, exponential, gamma, beta) can be written in a common form $f (x ∣ θ)$ $= h (x) exp (η (θ) T (x)$ $- A (θ)),$ for suitable functions $h$ , $η$ , $T$ and $A$ . For the Poisson, for example, $\frac{e^{- λ} λ^{x}}{x!}$ $= \frac{1}{x!} exp (x ln λ$ $- λ)$ , with $T (x) = x$ . This shared form is why they have so many properties in common: simple sufficient statistics (Chapter 8), conjugate Bayesian priors like the beta–binomial pair, and the unified "generalized linear models" used throughout applied statistics and machine learning (logistic regression is built on the Bernoulli family).

Why it matters: recognizing the family a distribution belongs to tells you, at once, many of the tools that will work for it.

5.9 Exercises

Each exercise has a hint and an answer. Draw a picture of the density and shade the area you want before calculating.

Exercises

A density by hand. The time $X$ (in hours) that a lab takes to process an urgent sample has density $f (x)$ $= c (4$ $- x)$ for $0 \leq x \leq 4$ . (a) Find $c$ . (b) Find $Pr (X \leq 1)$ . (c) Find $E [X]$ . (d) Find the median.

The density is a triangle with base 4 and height $4 c$ . For (d), the CDF is $F (x)$ $= c (4 x$ $- x^{2} / 2)$ ; set it equal to $\frac{1}{2}$ and solve the quadratic.

(a) $\int_{0}^{4} c (4$ $- x) d x$ $= 8 c$ $= 1$ , so $c = \frac{1}{8}$ . (b) $F (1)$ $= \frac{1}{8} (4 - \frac{1}{2})$ $= \frac{7}{16}$ $\approx 0.4375$ . (c) $\int_{0}^{4} x \frac{4 - x}{8} d x$ $= \frac{1}{8} (32 - \frac{64}{3})$ $= \frac{4}{3}$ hours. (d) $\frac{1}{8} (4 m$ $- m^{2} / 2)$ $= \frac{1}{2}$ gives $m^{2}$ $- 8 m$ $+ 8$ $= 0$ , so $m$ $= 4$ $- 2 \sqrt{2}$ $\approx 1.17$ hours.
Uniform. A random number generator produces $U$ $\sim Unif (0,$ $10)$ . (a) Find $Pr (3 \leq U \leq 7)$ . (b) Find $E [U]$ and $SD [U]$ . (c) Find $Pr (U = 5)$ , and explain why it isn't a problem that it equals what it does.

For a uniform distribution, probabilities are lengths divided by the total length.

(a) $\frac{7 - 3}{10} = 0.4$ . (b) Mean 5, SD $\sqrt{100 / 12}$ $\approx 2.89$ . (c) 0: a single point has no length. Every particular value has probability zero, but intervals don't, and some value always occurs.
Galton's daughters. Galton's 453 daughters had mean height 64.1 inches and SD 2.4 inches. Using a Gaussian model, find (a) the fraction taller than 68 inches, (b) the fraction between 62 and 66 inches, and (c) the height that only 10% of daughters exceed. (d) Check (a) against the data.

Standardize each value. For (d), with galton loaded as in Example 5.1, compute mean(galton$childHeight[galton$gender == "female"] >= 68).

(a) $z$ $= (68$ $- 64.1) / 2.4$ $= 1.625$ ; $1$ $- Φ (1.625)$ $\approx 0.052$ . (b) $Φ (0.792)$ $- Φ (- 0.875)$ $\approx 0.786$ $- 0.191$ $= 0.595$ . (c) $64.1$ $+ 1.2816$ $\times 2.4$ $\approx 67.2$ inches. (d) The data give 0.071, more than the model predicts, but almost all of the excess is heaping: 13 daughters were recorded at exactly 68 inches, which " $\geq 68$ " includes; only 4.2% were recorded above 68.
Low birth weight. In MASS::birthwt, birth weights have mean 2,945 g and SD 729 g. (a) Using a Gaussian model, estimate the fraction of babies under 2,500 g. (b) Compare with the data. (c) Why might a Gaussian model underestimate this fraction?

Standardize 2,500. For (b), mean(MASS::birthwt$bwt < 2500). For (c), think about the shape of the left tail: what kinds of babies are very light?

(a) $Φ (\frac{2500 - 2945}{729})$ $= Φ (- 0.61)$ $\approx 0.271$ . (b) The data give 0.312 (59 of 189). (c) Low weights come partly from a separate group (premature births) added to the main bell, so the left tail is heavier than a single Gaussian allows. Also, this hospital sample over-represented high-risk pregnancies. Models are often weakest in the tails, where the medically important cases are.
IQ. IQ scores are designed to be $N (100, 15^{2})$ . (a) What fraction of people score above 130? (b) What score puts someone in the top 1%? (c) In a city of 8 million, about how many people would you expect to score above 160?

(a) $z = 2$ ; (b) qnorm(0.99, 100, 15); (c) multiply 8 million by $Pr (Z > 4)$ .

(a) $1$ $- Φ (2)$ $\approx 0.023$ . (b) $100$ $+ 2.326$ $\times 15$ $\approx 134.9$ . (c) $Pr (Z > 4)$ $\approx 3.2$ $\times 10^{- 5}$ , so about $8,000,000$ $\times 0.000032$ $\approx 250$ people, though IQ tests aren't designed to measure accurately that far into the tail.
Exponential waits. Customers arrive at a pharmacy's prescription counter at an average rate of 12 per hour, independently. (a) What is the distribution of the time between arrivals, and its mean in minutes? (b) What is the probability of a gap of more than 10 minutes? (c) What is the median gap?

The gap is $Exp (12)$ in hours, or $Exp (0.2)$ in minutes.

(a) $Exp (0.2)$ per minute; mean 5 minutes. (b) $e^{- 0.2 \times 10}$ $= e^{- 2}$ $\approx 0.135$ . (c) $ln 2 / 0.2$ $\approx 3.47$ minutes.
Carbon dating. A wooden tool contains 30% of the carbon-14 that living wood contains. Using a half-life of 5,730 years, how old is it?

Solve ${(\frac{1}{2})}^{t / 5730}$ $= 0.30$ by taking logs.

$t$ $= 5730$ $\times \frac{ln 0.30}{ln 0.5}$ $\approx 5730$ $\times 1.737$ $\approx 9,950$ years.
Memorylessness. A website's server crashes at random, at an average rate of once every 30 days, as a Poisson process. It has now run for 45 days without a crash. (a) What is the probability it runs at least another 30 days? (b) A manager says: "It's been 45 days, so a crash is overdue." Explain the mistake. (c) When would the manager's reasoning be correct?

The time to a crash is exponential with mean 30 days. Use the memoryless property.

(a) $e^{- 30 / 30}$ $= e^{- 1}$ $\approx 0.368$ , the same as for a brand-new server. (b) Under the exponential model, the past 45 days are irrelevant; nothing is "due" (the gambler's fallacy again). (c) If crashes came from wear (memory leaks, filling disks), the failure rate would rise with time since the last restart, and a Weibull with shape $k > 1$ would be a better model; then a long run would make a crash more likely.
Gamma. Using Example 5.11's setting (911 calls at 4.5 per hour), find (a) the expected waiting time until the 5th call, and (b) the probability that the 5th call arrives within half an hour, in two ways.

$T$ $\sim Gamma (5,$ $4.5)$ . Alternatively, the 5th call arrives within half an hour if there are at least 5 calls in that half hour, a $Pois (2.25)$ count.

(a) $5 / 4.5 \approx 1.11$ hours, about 67 minutes. (b) pgamma(0.5, 5, 4.5) $=$ 1 - ppois(4, 2.25) $\approx 0.078$ .
Beta update. A forensic lab has tested 20 randomly chosen pills from a seizure and found that 12 contain fentanyl. (a) Starting from a uniform prior for the proportion $p$ of fentanyl pills, what is the posterior distribution? (b) Find the posterior mean. (c) Find the posterior probability that more than 60% of the pills contain fentanyl.

Beta prior $(1, 1)$ plus 12 successes and 8 failures. Use 1 - pbeta(0.6, a, b).

(a) $Beta (13, 9)$ . (b) $\frac{13}{22} \approx 0.59$ . (c) 1 - pbeta(0.6, 13, 9) $\approx 0.48$ : given the data, it's roughly a coin flip whether more than 60% of the pills contain fentanyl. Twenty pills can't pin the proportion down.
Log-Gaussian. Blood concentrations of a drug, one hour after a standard dose, are log-Gaussian with median 4 mg/L, and their natural logarithms have SD 0.5. (a) What is $μ$ for the log concentrations? (b) What fraction of patients exceed 10 mg/L, a toxic level? (c) Why is the mean concentration higher than the median?

The median of a log-Gaussian is $e^{μ}$ . Then $Pr (X > 10)$ $= 1$ $- Φ (\frac{ln 10 - μ}{0.5})$ .

(a) $μ$ $= ln 4$ $\approx 1.386$ . (b) $z$ $= (2.303$ $- 1.386) / 0.5$ $\approx 1.83$ ; $1$ $- Φ (1.83)$ $\approx 0.034$ , about 3.4% of patients. (c) The distribution is right-skewed: a few patients with very high concentrations pull the mean ( $e^{μ + σ^{2} / 2}$ $= 4 e^{0.125}$ $\approx 4.53$ ) above the median.
Simulation by inverse transform. (a) Show that if $U \sim Unif (0, 1)$ , then $X$ $= - \frac{1}{λ} ln (1$ $- U)$ is $Exp (λ)$ . (b) Use this in R, without rexp(), to simulate 10,000 waits with $λ = 2$ , and check their mean.

Solve $u$ $= F (x)$ $= 1$ $- e^{- λ x}$ for $x$ . In R: u <- runif(10000); x <- -log(1 - u) / 2.

(a) $u = 1 - e^{- λ x}$ gives $x$ $= - \frac{1}{λ} ln (1$ $- u)$ , which is $F^{- 1} (u)$ ; by the inverse transform theorem, $F^{- 1} (U)$ has CDF $F$ . (b) mean(-log(1 - runif(10000)) / 2) is close to $1 / 2$ .
Gaussian approximation. A coin is flipped 400 times. (a) Use the Gaussian approximation with continuity correction to estimate the probability of 220 or more heads. (b) Compare with the exact binomial. (c) Would 220 heads make you doubt that the coin is fair?

Mean 200, SD $\sqrt{400 \times 0.25}$ $= 10$ . "220 or more" becomes "above 219.5".

(a) $z$ $= (219.5$ $- 200) / 10$ $= 1.95$ ; $1$ $- Φ (1.95)$ $\approx 0.026$ . (b) 1 - pbinom(219, 400, 0.5) $\approx 0.026$ . (c) Somewhat: a fair coin gives 220 or more heads only about 2.6% of the time (or 5% for a result that far from 200 in either direction). Chapter 9 turns this reasoning into a formal test.

Chapter 5Continuous random variables

Think About It: Before You Start

Lecture 12 — Densities, expected value and the uniform

5.1 From histograms to densities

Galton's heights

Example 5.1: Galton's Sons ★ lecture example

Probability density functions

Definition: Continuous Random Variable and Density

Example 5.2: Where Does the Dart Land? ★ lecture example

The CDF, again

Example 5.3: The Median Dart ★ lecture example

Check Yourself

Think About It: How Precise Is a Measurement?

Key Takeaways

5.2 Expected value, variance and the uniform distribution

Definition: Expected Value and Variance of a Continuous Random Variable

Example 5.4: The Average Dart ★ lecture example

The uniform distribution

Definition: Uniform Distribution

Example 5.5: Waiting for the Bus ★ lecture example

Simulation: turning uniform numbers into anything

Think About It: The Inspection Paradox in the Wild

Key Takeaways

Lecture 13 — The Gaussian distribution

5.3 The Gaussian distribution

The bell curve

Definition: Gaussian Distribution

History: Why Gaussian and Not Normal?

Standardizing

Theorem: Standardizing a Gaussian Random Variable

Idea: The 68–95–99.7 Rule

Example 5.6: How Tall Were Galton's Sons? ★ lecture example

History: Galton's Heights, and Why He Collected Them

Causal Inference: Why Are So Many Things Gaussian?

Example 5.7: Just Over the Limit? ★ lecture example

The Gaussian approximation to the binomial

Example 5.8: The Airline, One Last Time

History: Ceres, Gauss and the Birth of the Bell Curve

Check Yourself

Think About It: Abnormal, or Just Unlucky?

Key Takeaways

Lecture 14 — Exponential, gamma, beta and other shapes

5.4 Waiting times: exponential and gamma

The exponential distribution

Definition: Exponential Distribution

Theorem: Waiting Times in a Poisson Process

Theorem: The Exponential Is Memoryless

Example 5.9: British Coal-Mining Disasters ★ lecture example

Example 5.10: How Old Is the Iceman? ★ lecture example

The gamma distribution

Definition: Gamma Distribution

Example 5.11: The Third Call

Check Yourself

Think About It: Memoryless or Not?

Key Takeaways

5.5 The beta distribution: uncertain proportions

Definition: Beta Distribution

Bayesian Inference: Learning a Free-Throw Percentage

Think About It: What Does the Prior Do?

5.6 More shapes of chance

Log-Gaussian: when causes multiply

Definition: log-Gaussian Distribution

Example 5.12: How Long Is a River? ★ lecture example

Weibull: things that wear out

Chi-squared and t: distributions for later

History: Student, Who Worked at a Brewery

Think About It: Which Scale Should We Predict on?

5.7 Choosing a continuous model

Idea: Algorithm: Which Continuous Model?

Check Yourself

Prediction and inference in this chapter

5.8 Interludes

Interlude (Optional): The Cauchy Distribution: When Averages Fail

Interlude (Optional): Families of Distributions

5.9 Exercises

Exercises

Chapter 5
Continuous random variables