Chapter 4
Discrete random variables

List the outcomes. Write R for a right answer and W for a wrong one. There are $2^3=8$ equally likely sequences (counting case 1 from Chapter 3), and $X$ counts the R's:

Group the outcomes by the value of $X$. $X=0$ for one outcome, $X=1$ for three, $X=2$ for three, $X=3$ for one. Each outcome has probability $\frac{1}{8}$, so

As a formula. The number of sequences with exactly $x$ right answers is the number of ways to choose which $x$ of the 3 questions you got right, $\left(\genfrac{}{}{0ex}{}{3}{x}\right)$ (counting case 3). So $p(x)$$=\left(\genfrac{}{}{0ex}{}{3}{x}\right){\left(\frac{1}{2}\right)}^3,$$$$x$$=0,$$1,$$2,$$3.$ (Check: $\left(\genfrac{}{}{0ex}{}{3}{1}\right)/8=3/8$.) This is our first binomial distribution; Section 4.4 generalizes it.

Valid? Every $p(x)\ge 0$, and $\frac{1}{8}$$+\frac{3}{8}$$+\frac{3}{8}$$+\frac{1}{8}$$=1$. Yes.

As a picture, and a check in R (dbinom computes binomial probabilities, as we'll see in Section 4.4):

x <- 0:3
p <- dbinom(x, size = 3, prob = 0.5)
p
sum(p)
ggplot(data.frame(x, p), aes(x, p)) +
  geom_col(width = 0.4) +
  labs(x = "number of correct guesses, x", y = "p(x)")

[1] 0.125 0.375 0.375 0.125
[1] 1

The bars for a PMF are drawn thin and separated, to show that $X$ can only take these exact values. Notice the good news and the bad news: a pure guesser gets at least 2 out of 3 half the time, but all 3 only one time in 8. (A teacher who wants to tell guessers from students who know the material needs more than three questions: see Exercise 13.)

The empirical PMF is the table of relative frequencies (Section 3.1): for each value $x$, the fraction of the 100 years with exactly $x$ discoveries.

d <- as.numeric(discoveries)       # 100 yearly counts, 1860 to 1959
table(d)
p_hat <- prop.table(table(d))       # relative frequencies
round(p_hat, 2)

d
 0  1  2  3  4  5  6  7  8  9 10 12 
 9 12 26 20 12  7  6  4  1  1  1  1 
d
   0    1    2    3    4    5    6    7    8    9   10   12 
0.09 0.12 0.26 0.20 0.12 0.07 0.06 0.04 0.01 0.01 0.01 0.01 

With exactly 100 years, each count divided by 100 is its estimated probability: $\hat{p}(2)=0.26$, the most common value, and so on. The counts add to 100, so the estimates add to 1.

ggplot(data.frame(x = as.numeric(names(p_hat)), p = as.numeric(p_hat)), aes(x, p)) +
  geom_col(width = 0.4) +
  scale_x_continuous(breaks = 0:12) +
  labs(x = "great discoveries in a year, x", y = "estimated p(x)")

The probability. $\{X\ge 4\}$ is the disjoint union of $\{X=4\},$$\{X=5\},$$\dots$, so add the estimates: $\widehat{\mathrm{\text{Pr}}}(X\ge 4)$$=\frac{12+7+6+4+1+1+1+1}{100}$$=\frac{33}{100}$$=0.33.$ In R: mean(d >= 4).

A critical reading. What counts as "great" is someone's judgement, made decades later, and recent discoveries may not yet have been recognized: the measurement process (Section 1.4) shapes the data. Notice also $\hat{p}(11)=0$: no year had exactly 11, though 12 happened once. An empirical PMF with zeros in the middle reminds us that it's an estimate of the true distribution, not the distribution itself.

Running totals. Add the PMF from left to right:

Between the possible values, $F$ stays flat. So in full, $$F(x)=\{\begin{array}{cc}0,\hfill & x<0,\hfill \\ \frac{1}{8},\hfill & 0\le x<1,\hfill \\ \frac{1}{2},\hfill & 1\le x<2,\hfill \\ \frac{7}{8},\hfill & 2\le x<3,\hfill \\ 1,\hfill & x\ge 3.\hfill \end{array}$$ Notice that $F(1.5)=\frac{1}{2}$, for example: "at most 1.5 right" means "at most 1 right".

ggplot(data.frame(x = c(-1, 0:3, 4)), aes(x)) +
  stat_function(fun = function(x) pbinom(x, 3, 0.5), geom = "step", n = 1001) +
  geom_point(data = data.frame(x = 0:3, F = pbinom(0:3, 3, 0.5)), aes(x, F)) +
  labs(x = "x", y = "F(x) = Pr(X ≤ x)")

The probabilities.

• $\mathrm{\text{Pr}}(X\le 1)$$=F(1)$$=\frac{1}{2}$.
• $\mathrm{\text{Pr}}(X<1)$$=\mathrm{\text{Pr}}(X\le 0)$$=F(0)$$=\frac{1}{8}$.
• $\mathrm{\text{Pr}}(X\ge 2)$$=1$$-\mathrm{\text{Pr}}(X\le 1)$$=1$$-F(1)$$=\frac{1}{2}$. So a pure guesser passes half the time.
• $\mathrm{\text{Pr}}(1<X\le 3)$$=F(3)$$-F(1)$$=1$$-\frac{1}{2}$$=\frac{1}{2}$. (Check: that's $\mathrm{\text{Pr}}(X=2)$$+\mathrm{\text{Pr}}(X=3)$$=\frac{3}{8}$$+\frac{1}{8}$.)

In R, pbinom(1, 3, 0.5) gives $F(1)$. Every named distribution has a "p" function for its CDF.

By hand. Multiply each value by its estimated probability and add: $$\begin{array}{cc}\hfill \mathrm{\text{E}}[X]& =0(0.09)+1(0.12)+2(0.26)+3(0.20)+4(0.12)+5(0.07)\hfill \\ \hfill & +6(0.06)+7(0.04)+8(0.01)+9(0.01)+10(0.01)+12(0.01)\hfill \\ \hfill & =0+0.12+0.52+0.60+0.48+0.35+0.36+0.28+0.08+0.09+0.10+0.12\hfill \\ \hfill & =3.10.\hfill \end{array}$$$\mathrm{\text{E}}[X]$$=0(0.09)$$+1(0.12)$$+2(0.26)$$+3(0.20)$$+4(0.12)$$+5(0.07)$$$$+6(0.06)$$+7(0.04)$$+8(0.01)$$+9(0.01)$$+10(0.01)$$+12(0.01)$$=0$$+0.12$$+0.52$$+0.60$$+0.48$$+0.35$$+0.36$$+0.28$$+0.08$$+0.09$$+0.10$$+0.12$$=3.10.$

In R, the same calculation is sum(x * p). Compare it with the ordinary average of the data:

x <- as.numeric(names(p_hat))
sum(x * p_hat)          # expected value of the empirical PMF
mean(d)                 # the average of the 100 yearly counts

[1] 3.1
[1] 3.1

Both give 3.1. That's no coincidence: averaging 100 numbers is the same as weighting each distinct value by the fraction of the time it occurs, which is exactly what $\sum x\hat{p}(x)$ does. The expected value of an empirical distribution is the sample mean of Chapter 2. For a model distribution, $\mathrm{\text{E}}[X]$ is the mean we'd see in the long run.

Meaning. Averaged over the century, about 3.1 great discoveries were made per year: 310 in 100 years. No year had 3.1 discoveries, but a historian estimating how many great discoveries to expect in a decade would say about 31. Try the "great discoveries" option in the widget above: the running average of simulated years settles at 3.1.

Prediction and inference. As a prediction, 3.1 is our best guess (under squared loss) for the number of great discoveries in a year like these. As an inference, it estimates the underlying rate of discovery in that century. Notice that the two uses need different caveats: the prediction assumes future years behave like 1860–1959; the inference assumes the years are a fair record of what happened.

Set up. Let $W$ be the prize on one ticket. Each tier's probability is 1 divided by its odds, and the tiers are disjoint (a ticket wins at most one), so by the definition of expected value $$\mathrm{\text{E}}[W]=\sum _{\mathrm{\text{tiers}}}(\mathrm{\text{prize}})\times \mathrm{\text{Pr}}(\mathrm{\text{tier}}).$$

(a) The non-jackpot tiers contribute, for example, $4/38.32$$\approx \$0.104$ for the first and $\mathrm{1,000,000}/\mathrm{11,688,053.52}$$\approx \$0.086$ for the $1 million tier. Adding all eight in R:

prize <- c(4, 4, 7, 7, 100, 100, 50000, 1e6)
odds  <- c(38.32, 91.98, 701.33, 579.76, 14494.11, 36525.17, 913129.18, 11688053.52)
round(prize / odds, 4)         # each tier's contribution, in dollars
sum(prize / odds)              # expected non-jackpot prize per ticket
(2 - sum(prize / odds)) * 292201338    # break-even jackpot

[1] 0.1044 0.0435 0.0100 0.0121 0.0069 0.0027 0.0548 0.0856
[1] 0.3198782
[1] 490933824

So the non-jackpot prizes are worth about 32 cents per $2 ticket. Interestingly, the $50,000 and $1 million tiers contribute 14 cents of that between them: the rare big prizes matter about as much as the frequent small ones.

(b) The jackpot adds $J/\mathrm{292,201,338}$. For $\mathrm{\text{E}}[W]=2$ we need $0.32$$+\frac{J}{\mathrm{292,201,338}}$$=2$$\Rightarrow$$J$$\approx 1.68$$\times \mathrm{292,201,338}$$\approx \$491\mathrm{\text{million}}.$

(c) Three reasons a $491 million advertised jackpot still isn't a fair bet. The advertised jackpot is paid over 29 years; the lump-sum cash value is roughly half. Large prizes are taxed heavily. And when the jackpot is huge, ticket sales soar, so the chance of having to share the jackpot with other winners rises. Each of these cuts the expected value well below $2. On expected value alone, a Powerball ticket is a poor investment; people buy them for the dream, which is a perfectly reasonable thing to pay $2 for, as long as you know that's what you're buying.

Source: Powerball prize chart, powerball.com.

Where the expected value comes from. A bar chart of each tier's contribution (prize × probability) shows that the rare large prizes contribute as much as the common small ones:

tiers <- c("PB only", "1+PB", "2+PB", "3", "3+PB", "4", "4+PB", "5")
d <- data.frame(tier = factor(tiers, levels = tiers), cents = 100 * prize / odds)
ggplot(d, aes(tier, cents)) + geom_col(width = 0.6) +
  labs(x = "prize tier (numbers matched)", y = "cents per ticket")

Set up. Let $X$ be the number who show up; its PMF is binomial (Section 4.4), and R's dbinom(x, 110, 0.9) gives $\mathrm{\text{Pr}}(X=x)$. The cost is a function of $X$: $c(x)$$=\{\begin{array}{cc}250(100-x),\hfill & x<100,\hfill \\ 1000(x-100),\hfill & x\ge 100.\hfill \end{array}$

Apply the theorem. $\mathrm{\text{E}}[c(X)]$$={\sum }_{x=0}^{110}c(x)\mathrm{\text{Pr}}(X=x)$, a sum of 111 terms that R does in one line:

x <- 0:110
p <- dbinom(x, size = 110, prob = 0.9)
cost <- ifelse(x < 100, 250 * (100 - x), 1000 * (x - 100))
sum(cost * p)                       # expected cost per flight
# repeat for every number of tickets from 95 to 120
exp_cost <- sapply(95:120, function(n) {
  x <- 0:n
  sum(ifelse(x < 100, 250 * (100 - x), 1000 * (x - 100)) * dbinom(x, n, 0.9))
})
(95:120)[which.min(exp_cost)]       # the best number of tickets

[1] 1240.143
[1] 108

The exact expected cost is $1240.14, matching the simulation's estimate. Doing the same for every number of tickets shows that 108 tickets minimizes the expected cost (about $1041 per flight), as the widget in Chapter 1 showed.

ggplot(data.frame(n = 95:120, exp_cost), aes(n, exp_cost)) +
  geom_line() + geom_point(size = 1) +
  geom_point(data = data.frame(n = 108, exp_cost = min(exp_cost)), size = 3.5, colour = "orange") +
  labs(x = "tickets sold", y = "expected cost ($)")

Chapter 1's question ("how many tickets?") is now a precise calculation: choose the decision that minimizes the expected value of the loss.

Expected values. $$\begin{array}{cc}\hfill \mathrm{\text{E}}[A]& =\mathrm{50,000},\hfill \\ \hfill \mathrm{\text{E}}[B]& =\mathrm{100,000}\times \frac{1}{2}+0\times \frac{1}{2}=\mathrm{50,000},\hfill \\ \hfill \mathrm{\text{E}}[C]& =\mathrm{200,000}\times \frac{1}{4}+0\times \frac{3}{4}=\mathrm{50,000}.\hfill \end{array}$$$\mathrm{\text{E}}[A]$$=\mathrm{50,000},$$\mathrm{\text{E}}[B]$$=\mathrm{100,000}$$\times \frac{1}{2}$$+0$$\times \frac{1}{2}$$=\mathrm{50,000},$$\mathrm{\text{E}}[C]$$=\mathrm{200,000}$$\times \frac{1}{4}$$+0$$\times \frac{3}{4}$$=\mathrm{50,000}.$ All three have the same expected value. On average, over many plays, they would pay the same.

Variances, by the shortcut formula $\mathrm{\text{Var}}$$=\mathrm{\text{E}}[X^2]$$-(\mathrm{\text{E}}[X]{)}^2$. Work in thousands of dollars to keep the numbers small:

• A always pays 50, so $\mathrm{\text{E}}[A^2]$$={50}^2$$=2500$ and $\mathrm{\text{Var}}[A]$$=2500$$-2500$$=0$: no risk at all.
• B: $\mathrm{\text{E}}[B^2]$$={100}^2$$\times \frac{1}{2}$$=5000$, so $\mathrm{\text{Var}}[B]$$=5000$$-2500$$=2500$ and $\mathrm{\text{SD}}[B]=50$.
• C: $\mathrm{\text{E}}[C^2]$$={200}^2$$\times \frac{1}{4}$$=\mathrm{10,000}$, so $\mathrm{\text{Var}}[C]$$=\mathrm{10,000}$$-2500$$=7500$ and $\mathrm{\text{SD}}[C]\approx 86.6$.

So the standard deviations are $0, $50,000 and about $86,600: the same average, increasingly risky.

sd_of <- function(x, p) { mu <- sum(x * p); sqrt(sum(x^2 * p) - mu^2) }
c(A = sd_of(50000, 1),
  B = sd_of(c(0, 100000), c(1/2, 1/2)),
  C = sd_of(c(0, 200000), c(3/4, 1/4)))

       A        B        C 
    0.00 50000.00 86602.54 

Which to choose? Expected value can't decide this: they're tied. Most people take A, because $50,000 for certain changes their life and a 75% chance of nothing (choice C) is painful. A billionaire, for whom $50,000 is small change, might reasonably take C for the fun of it. The standard deviation measures the risk; how much risk you should take depends on your circumstances, which is the loss function idea of Chapter 2 again. The real game show Deal or No Deal offered contestants exactly this kind of choice, and economists studied the contestants' decisions to learn how people value risk.

The three choices side by side. Same balance point, very different spread:

choices <- data.frame(
  choice = rep(c("A", "B", "C"), c(1, 2, 2)),
  prize = c(50, 0, 100, 0, 200), prob = c(1, 0.5, 0.5, 0.75, 0.25))
ggplot(choices, aes(prize, prob)) + geom_col(width = 12) +
  geom_vline(xintercept = 50, linetype = "dashed") +
  facet_wrap(~ choice) + scale_x_continuous(breaks = c(0, 100, 200)) +
  theme(panel.spacing = unit(1.2, "lines")) +
  labs(x = "prize (thousands of dollars)", y = "probability")

Is it binomial? Each prediction is right or wrong (binary); a random chooser's guesses are independent; the number of matches is fixed; and each guess is right with probability $\frac{1}{2}$ (same $p$). So the number right is binomial.

(a) $X$$\sim \mathrm{\text{Bin}}(8,\frac{1}{2})$: $\mathrm{\text{Pr}}(X=8)$$={\left(\frac{1}{2}\right)}^8$$=\frac{1}{256}$$\approx 0.0039.$

(b) $Y$$\sim \mathrm{\text{Bin}}(14,\frac{1}{2})$. "12 or more" is $Y=12,13$ or $14$: $\mathrm{\text{Pr}}(Y\ge 12)$$=\frac{\left(\genfrac{}{}{0ex}{}{14}{12}\right)+\left(\genfrac{}{}{0ex}{}{14}{13}\right)+\left(\genfrac{}{}{0ex}{}{14}{14}\right)}{2^{14}}$$=\frac{91+14+1}{\mathrm{16,384}}$$=\frac{106}{\mathrm{16,384}}$$\approx 0.0065.$

dbinom(8, 8, 0.5)                    # (a)
1 - pbinom(11, 14, 0.5)              # (b) Pr(Y >= 12)
200 * dbinom(8, 8, 0.5)              # (c) expected number of perfect animals

[1] 0.00390625
[1] 0.006469727
[1] 0.78125

(c) Each animal is perfect with probability $\frac{1}{256}$, so by linearity (with indicators, as in Example 4.7) the expected number of perfect animals is $200$$\times \frac{1}{256}$$\approx 0.78$. With a few hundred animals predicting, it's quite likely that some animal gets them all right by luck alone, and that animal is the one that becomes famous.

The lesson is the birthday problem again (Example 3.9): a one-in-256 event is rare for one animal and unremarkable among hundreds. We only heard about Paul because he got them right; nobody reported on the parakeet that got 4 out of 8. Selecting the winners after the fact is how many "amazing" track records, from stock-pickers to psychics, are made. Chapter 9 calls this the problem of multiple testing.

How surprising is 8 out of 8? The PMF of $\mathrm{\text{Bin}}(8,\frac{1}{2})$ shows how far Paul's result sits from what guessing usually produces:

d <- data.frame(x = 0:8, p = dbinom(0:8, 8, 0.5), paul = 0:8 == 8)
ggplot(d, aes(x, p, fill = paul)) + geom_col(width = 0.5) +
  scale_fill_manual(values = c("grey60", "orange"), guide = "none") +
  scale_x_continuous(breaks = 0:8) +
  labs(x = "correct predictions out of 8, if guessing", y = "probability")

This is inference. We asked: how probable is the observed result if the mechanism is pure guessing? Very improbable for one animal, but not for the luckiest of hundreds. Chapter 9 formalizes this as a hypothesis test and its p-value.

Is it binomial? Binary (positive or not); independent (stated); $n=20$ fixed; same $p=0.05$ for every clean athlete. So $X$ = number of false positives $\sim \mathrm{\text{Bin}}(20,0.05)$, with "success" meaning a positive test. (A "success" is just the outcome we count, here an unwelcome one.)

(a) $\mathrm{\text{Pr}}(X=0)$$=\left(\genfrac{}{}{0ex}{}{20}{0}\right){(0.05)}^0{(0.95)}^{20}$$={0.95}^{20}$$\approx 0.358.$

(b) $\mathrm{\text{Pr}}(X=1)$$=\left(\genfrac{}{}{0ex}{}{20}{1}\right){(0.05)}^1{(0.95)}^{19}$$=20$$\times 0.05$$\times {0.95}^{19}$$\approx 0.377.$

(c) "At least one" is easiest through the complement: $\mathrm{\text{Pr}}(X\ge 1)$$=1$$-\mathrm{\text{Pr}}(X=0)$$\approx 1$$-0.358$$=0.642.$

(d) $\mathrm{\text{E}}[X]$$=np$$=20$$\times 0.05$$=1$ and $\mathrm{\text{SD}}[X]$$=\sqrt{np(1-p)}$$=\sqrt{20\times 0.05\times 0.95}$$=\sqrt{0.95}$$\approx 0.97.$

dbinom(0, size = 20, prob = 0.05)       # (a)
dbinom(1, size = 20, prob = 0.05)       # (b)
1 - pbinom(0, size = 20, prob = 0.05)   # (c): 1 - Pr(X <= 0)
round(dbinom(0:5, 20, 0.05), 4)         # the first few PMF values

[1] 0.3584859
[1] 0.3773536
[1] 0.6415141
[1] 0.3585 0.3774 0.1887 0.0596 0.0133 0.0022

ggplot(data.frame(x = 0:6, p = dbinom(0:6, 20, 0.05)), aes(x, p)) +
  geom_col(width = 0.4) +
  labs(x = "false positives among 20 clean athletes", y = "probability")

What it means. Even with no drug users at all, there's nearly a two-in-three chance that someone tests positive. A single positive on a large team is not, by itself, evidence of drug use; that's why positive screens are confirmed with a second, more specific test (and why Section 3.6's Bayes calculation matters).