19 Random Variables

Chapter 18 assigned probabilities to events. But you rarely want the probability of an abstract subset — you want to know how tall someone is, how many clicks a page gets, how large the error will be.

A random variable attaches a number to each outcome, and once you have numbers you can average them, measure their spread, and add them up. That is what makes the rest of statistics possible.

19.1 What a random variable is

A random variable is a function from the sample space to the real numbers:

\[ X: \Omega \to \mathbb{R} \]

It is neither random nor a variable — it is a deterministic function (Section 3.4). The randomness is in which outcome occurs; \(X\) just reports a number once it has.

Experiment	A random variable on it
roll a die	\(X\) = the value shown
flip two coins	\(X\) = number of heads
pick a person	\(X\) = their height
fit a model	\(X\) = the prediction error

Convention: uppercase for the random variable, lowercase for a value it takes. So \(P(X = x)\) asks for the probability that the variable \(X\) takes the particular value \(x\).

19.2 Discrete and continuous

Discrete variables take values you can list — counts, categories, die faces — possibly infinitely many, as long as they are separated.

Continuous variables take any value in an interval. Heights, times, temperatures.

The distinction matters because probability behaves differently. For a continuous variable, \(P(X = x) = 0\) for every \(x\): there are uncountably many possible values, so no single one can carry positive probability. Only intervals do.

That sounds paradoxical — the variable does take some value — but it is the same fact as a point having zero length while an interval has positive length.

19.3 Probability mass and density functions

For a discrete variable, the probability mass function gives each value’s probability:

\[ p(x) = P(X = x), \qquad \sum_x p(x) = 1 \]

# X = number of heads in two fair coin flips
xs <- 0:2
pmf <- c(0.25, 0.5, 0.25)
rbind(x = xs, p = pmf, cumulative = cumsum(pmf))

           [,1] [,2] [,3]
x          0.00 1.00 2.00
p          0.25 0.50 0.25
cumulative 0.25 0.75 1.00

For a continuous variable, the probability density function gives probability per unit length, and you integrate it (Section 14.11):

\[ P(a \leq X \leq b) = \int_a^b f(x)\,dx, \qquad \int_{-\infty}^{\infty} f(x)\,dx = 1 \]

Watch out

A density is not a probability. It can exceed 1 — a uniform distribution on \([0, 0.5]\) has density 2 everywhere on that interval. What must be at most 1 is the integral, not the height.

Only after integrating over an interval do you get something interpretable as a probability.

# density of Uniform(0, 0.5) is 2, yet the total is 1
c(density_height = 1 / 0.5, total_area = 2 * 0.5)

density_height     total_area 
             2              1

19.4 The cumulative distribution function

The CDF works for both kinds and is often the more convenient object:

\[ F(x) = P(X \leq x) \tag{19.1}\]

Every CDF is non-decreasing, approaches 0 at \(-\infty\) and 1 at \(+\infty\). For discrete variables it is a staircase; for continuous ones a smooth climb, and there \(F'(x) = f(x)\) — the CDF and density are related by the fundamental theorem of calculus (Section 14.4).

draw_bar(
  x = c("0 heads", "1 head", "2 heads"),
  y = pmf,
  ylab = "probability"
)

Figure 19.1: The probability mass function for the number of heads in two coin flips. One head is twice as likely as either extreme, because there are two ways to get it.

xg <- seq(-0.5, 2.5, length.out = 400)
Fx <- sapply(xg, function(v) sum(pmf[xs <= v]))
draw_line(
  x = xg,
  y = Fx,
  points = FALSE,
  xlab = "x",
  ylab = "F(x) = P(X <= x)"
)

Figure 19.2: The matching CDF. It is a staircase: flat between attainable values, jumping by \(p(x)\) at each one, and the jump heights are exactly the mass function.

Probabilities of intervals come straight off the CDF, which is why it is what software computes — pnorm(), pbinom() and friends are all CDFs:

\[ P(a < X \leq b) = F(b) - F(a) \]

19.5 Expectation

The expectation is the probability-weighted average — the long-run mean:

\[ \mathbb{E}[X] = \sum_x x\,p(x) \qquad\text{or}\qquad \mathbb{E}[X] = \int_{-\infty}^{\infty} x\,f(x)\,dx \tag{19.2}\]

die <- 1:6
p_die <- rep(1 / 6, 6)
EX <- sum(die * p_die)
EX

[1] 3.5

\(3.5\) — a value the die can never show. The expectation need not be attainable; it is a summary of the distribution, not a prediction of any single outcome.

The property that makes expectation indispensable is linearity:

\[ \mathbb{E}[aX + b] = a\,\mathbb{E}[X] + b \qquad \mathbb{E}[X + Y] = \mathbb{E}[X] + \mathbb{E}[Y] \tag{19.3}\]

The second holds whether or not \(X\) and \(Y\) are independent, which is unusual and extremely useful — it follows from the linearity of sums and integrals (Section 14.6), nothing more.

To take the expectation of a function of \(X\), weight the transformed values:

\[ \mathbb{E}[g(X)] = \sum_x g(x)\,p(x) \]

c(
  E_X = EX,
  E_X_squared = sum(die^2 * p_die),
  E_2X_plus_1 = sum((2 * die + 1) * p_die),
  check_linearity = 2 * EX + 1
)

            E_X     E_X_squared     E_2X_plus_1 check_linearity 
        3.50000        15.16667         8.00000         8.00000

Watch out

\(\mathbb{E}[g(X)] \neq g(\mathbb{E}[X])\) in general. Here \(\mathbb{E}[X^2] = 15.17\) while \((\mathbb{E}[X])^2 = 12.25\).

For a convex \(g\) the inequality always runs one way — Jensen’s inequality says \(\mathbb{E}[g(X)] \geq g(\mathbb{E}[X])\) — which is Section 15.3 applied to a distribution. The gap here is \(2.92\), and it is exactly the variance.

19.6 Variance and standard deviation

The variance measures spread as the average squared distance from the mean:

\[ \operatorname{Var}(X) = \mathbb{E}\bigl[(X - \mu)^2\bigr] = \mathbb{E}[X^2] - \bigl(\mathbb{E}[X]\bigr)^2 \tag{19.4}\]

The second form is usually easier to compute and is the identity from Section 14.13.

VarX <- sum(die^2 * p_die) - EX^2
c(variance = VarX, as_fraction = 35 / 12, sd = sqrt(VarX))

   variance as_fraction          sd 
   2.916667    2.916667    1.707825

The standard deviation \(\sigma = \sqrt{\operatorname{Var}(X)}\) is preferred for reporting because it has the same units as \(X\); a variance of “2.92 squared pips” is hard to interpret.

Scaling behaves differently from the mean:

\[ \operatorname{Var}(aX + b) = a^2\operatorname{Var}(X) \]

Adding a constant shifts the distribution without changing its spread, so \(b\) drops out. Multiplying by \(a\) scales distances by \(a\) and squared distances by \(a^2\).

For sums, the covariance term appears unless the variables are uncorrelated:

\[ \operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y) + 2\operatorname{Cov}(X,Y) \tag{19.5}\]

19.7 Moments

The \(k\)-th moment is \(\mathbb{E}[X^k]\), and the \(k\)-th central moment is \(\mathbb{E}[(X-\mu)^k]\). The first few have names and describe shape:

Moment	Name	Describes
1st	mean	location
2nd central	variance	spread
3rd standardized	skewness	asymmetry
4th standardized	kurtosis	tail weight

Skewness is zero for a symmetric distribution, positive when the right tail is longer. Kurtosis measures how much probability sits in the tails, and heavy tails are what make extreme events more common than a normal model suggests.

set.seed(2)
z <- rnorm(20000)
e <- rexp(20000)
skew <- function(x) mean(((x - mean(x)) / sd(x))^3)
rbind(
  normal = c(skew = skew(z), sd = sd(z)),
  exponential = c(skew = skew(e), sd = sd(e))
)

                   skew       sd
normal      -0.01366767 1.005924
exponential  1.91180253 1.003259

The normal is symmetric so its skewness is near zero; the exponential has a long right tail and skewness near 2.

Watch out

Moments need not exist. If the defining integral diverges (Section 14.7), the moment is undefined. The Cauchy distribution has no mean at all — its sample average does not settle down no matter how much data you collect, and the law of large numbers simply does not apply.

Heavy-tailed data in the wild is not always that extreme, but it is common enough that “the mean of my sample” can be a meaningless summary.

19.8 Transformations of random variables

Applying a function to a random variable gives another random variable, and the distribution changes in ways that are easy to get wrong.

For a linear transformation everything is simple:

\[ Y = aX + b \;\Longrightarrow\; \mathbb{E}[Y] = a\mu + b, \quad \operatorname{Var}(Y) = a^2\sigma^2 \]

The most common case is standardization, which centers and rescales to mean 0 and variance 1:

\[ Z = \frac{X - \mu}{\sigma} \tag{19.6}\]

zdie <- (die - EX) / sqrt(VarX)
c(
  mean = sum(zdie * p_die),
  variance = sum(zdie^2 * p_die) - sum(zdie * p_die)^2
)

    mean variance 
       0        1

For non-linear transformations there is no such shortcut, and the naive guess is wrong: \(\mathbb{E}[1/X] \neq 1/\mathbb{E}[X]\), \(\mathbb{E}[\log X] \neq \log\mathbb{E}[X]\). Jensen’s inequality tells you the direction but not the size.

19.9 Joint distributions

Two variables together are described by a joint distribution \(p(x,y) = P(X=x, Y=y)\). Summing out one variable recovers the other’s marginal distribution:

\[ p_X(x) = \sum_y p(x, y) \]

joint <- matrix(
  c(0.1, 0.2, 0.3, 0.4),
  nrow = 2, byrow = TRUE,
  dimnames = list(X = c("0", "1"), Y = c("0", "1"))
)
joint

   Y
X     0   1
  0 0.1 0.2
  1 0.3 0.4

list(
  marginal_X = rowSums(joint),
  marginal_Y = colSums(joint)
)

$marginal_X
  0   1 
0.3 0.7 

$marginal_Y
  0   1 
0.4 0.6

The word marginal is literal: these are the row and column sums, traditionally written in the margins of the table.

19.10 Covariance and correlation

Covariance measures whether two variables move together:

\[ \operatorname{Cov}(X, Y) = \mathbb{E}\bigl[(X-\mu_X)(Y-\mu_Y)\bigr] = \mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y] \tag{19.7}\]

Positive when they tend to be large together, negative when one is large as the other is small. Note \(\operatorname{Cov}(X,X) = \operatorname{Var}(X)\) — variance is covariance with itself.

Its problem is units: covariance is measured in (units of \(X\))\(\times\)(units of \(Y\)), so its size means nothing on its own. Correlation fixes that by standardizing:

\[ \rho = \frac{\operatorname{Cov}(X,Y)}{\sigma_X\sigma_Y} \in [-1, 1] \tag{19.8}\]

This is exactly the cosine similarity of Section 4.15, applied to centered variables — which is why it is bounded by \(\pm1\), for the same reason a cosine is.

lvl <- c(0, 1)
EXY <- sum(outer(lvl, lvl) * joint)
EXj <- sum(lvl * rowSums(joint))
EYj <- sum(lvl * colSums(joint))
c(
  E_XY = EXY, E_X = EXj, E_Y = EYj,
  covariance = EXY - EXj * EYj
)

      E_XY        E_X        E_Y covariance 
      0.40       0.70       0.60      -0.02

set.seed(9)
n <- 3000
a <- rnorm(n)
b <- 0.6 * a + 0.8 * rnorm(n)
c(
  cov = cov(a, b),
  cor = cor(a, b),
  var_sum = var(a + b),
  by_formula = var(a) + var(b) + 2 * cov(a, b)
)

       cov        cor    var_sum by_formula 
 0.5741221  0.5925460  3.0861047  3.0861047

Equation 19.5 holds exactly.

Watch out

Zero correlation does not mean independent. Correlation detects only linear association.

Take \(X\) uniform on \(\{-1, 0, 1\}\) and \(Y = X^2\). Then \(Y\) is a deterministic function of \(X\) — as dependent as two variables can be — yet:

xv <- c(-1, 0, 1)
pv <- rep(1 / 3, 3)
yv <- xv^2
c(
  E_X = sum(xv * pv),
  E_Y = sum(yv * pv),
  E_XY = sum(xv * yv * pv),
  covariance = sum(xv * yv * pv) - sum(xv * pv) * sum(yv * pv)
)

       E_X        E_Y       E_XY covariance 
 0.0000000  0.6666667  0.0000000  0.0000000

Covariance exactly zero. The relationship is perfectly strong and perfectly non-linear, and correlation cannot see it.

xd <- runif(400, -2, 2)
yd <- xd^2 + rnorm(400, sd = 0.15)
draw_scatter(
  x = xd, y = yd,
  xlab = "x", ylab = "y"
)

Figure 19.3: A cloud with correlation near zero that is nonetheless completely determined: \(Y = X^2\) with noise. Any summary that reports only a correlation would call these variables unrelated.

round(cor(xd, yd), 4)

[1] -0.0163

19.11 Conditional expectation

The expectation of \(Y\) once you know \(X\):

\[ \mathbb{E}[Y \mid X = x] = \sum_y y\,p(y \mid x) \]

Crucially, \(\mathbb{E}[Y \mid X]\) is itself a random variable — a function of \(X\), which is random. Averaging it recovers the unconditional mean:

\[ \mathbb{E}\bigl[\mathbb{E}[Y \mid X]\bigr] = \mathbb{E}[Y] \tag{19.9}\]

the law of total expectation — Equation 18.5 with expectations in place of probabilities.

In machine learning

Conditional expectation is what regression estimates. The function minimizing expected squared error is exactly \(\mathbb{E}[Y \mid \mathbf{X} = \mathbf{x}]\), so a regression model is an estimate of a conditional expectation — which is why squared error is the default loss, and why a model fit that way predicts the mean rather than the median or the mode.

19.12 Independence of random variables

\(X\) and \(Y\) are independent when their joint factors into the product of marginals:

\[ p(x, y) = p_X(x)\,p_Y(y) \quad\text{for all } x, y \tag{19.10}\]

which extends Equation 18.4 from events to variables.

product_table <- outer(rowSums(joint), colSums(joint))
round(product_table, 4)

     0    1
0 0.12 0.18
1 0.28 0.42

all.equal(as.vector(joint), as.vector(product_table))

[1] "Mean relative difference: 0.08"

The joint does not factor, so these are dependent — consistent with the non-zero covariance computed earlier.

Independence implies more than zero covariance:

If \(X \perp Y\)	Statement
covariance	\(\operatorname{Cov}(X,Y) = 0\)
products	\(\mathbb{E}[XY] = \mathbb{E}[X]\mathbb{E}[Y]\)
variance of a sum	\(\operatorname{Var}(X+Y) = \operatorname{Var}(X)+\operatorname{Var}(Y)\)
any functions	\(g(X)\) and \(h(Y)\) are independent too

The implications run one way only. Independence gives zero covariance; zero covariance does not give independence, as the \(Y = X^2\) example showed.

19.13 Summary

Concept	Definition
Random variable	a function \(\Omega \to \mathbb{R}\)
PMF / PDF	\(P(X=x)\) / density, integrates to 1
CDF	\(F(x) = P(X \leq x)\)
Expectation	\(\mathbb{E}[X] = \sum x\,p(x)\)
Linearity	\(\mathbb{E}[X+Y] = \mathbb{E}[X]+\mathbb{E}[Y]\), always
Variance	\(\mathbb{E}[X^2] - (\mathbb{E}[X])^2\)
Scaling	\(\operatorname{Var}(aX+b) = a^2\operatorname{Var}(X)\)
Covariance	\(\mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]\)
Correlation	\(\operatorname{Cov}/(\sigma_X\sigma_Y) \in [-1,1]\)
Independence	\(p(x,y) = p_X(x)p_Y(y)\)

19.14 Exercises

1. For \(X\) = the number of heads in three fair coin flips, find the PMF, \(\mathbb{E}[X]\) and \(\operatorname{Var}(X)\).

Solution

x3 <- 0:3
p3 <- choose(3, x3) / 2^3
E3 <- sum(x3 * p3)
V3 <- sum(x3^2 * p3) - E3^2
rbind(x = x3, p = p3)

   [,1]  [,2]  [,3]  [,4]
x 0.000 1.000 2.000 3.000
p 0.125 0.375 0.375 0.125

c(mean = E3, variance = V3)

    mean variance 
    1.50     0.75

Mean \(1.5\) and variance \(0.75\). These match the binomial formulas \(np\) and \(np(1-p)\) with \(n=3\), \(p=0.5\) — which Section 20.2 derives.

2. A fair die is rolled. Let \(Y = 2X + 1\). Find \(\mathbb{E}[Y]\) and \(\operatorname{Var}(Y)\) two ways.

Solution

y <- 2 * die + 1
c(
  E_direct = sum(y * p_die),
  E_formula = 2 * EX + 1,
  Var_direct = sum(y^2 * p_die) - sum(y * p_die)^2,
  Var_formula = 4 * VarX
)

   E_direct   E_formula  Var_direct Var_formula 
    8.00000     8.00000    11.66667    11.66667

The variance quadruples while the mean only doubles-and-shifts — the \(a^2\) in \(\operatorname{Var}(aX+b)\), and the reason standard deviation (which scales by \(|a|\)) is the more intuitive summary.

3. Show that \(\operatorname{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2\) from the definition.

Solution

Expand the square and use linearity:

\[ \begin{aligned} \mathbb{E}[(X-\mu)^2] &= \mathbb{E}[X^2 - 2\mu X + \mu^2] \\ &= \mathbb{E}[X^2] - 2\mu\,\mathbb{E}[X] + \mu^2 \\ &= \mathbb{E}[X^2] - 2\mu^2 + \mu^2 = \mathbb{E}[X^2] - \mu^2 \end{aligned} \]

c(
  definition = sum((die - EX)^2 * p_die),
  shortcut = sum(die^2 * p_die) - EX^2
)

definition   shortcut 
  2.916667   2.916667

The step pulling \(\mu\) out of \(\mathbb{E}[2\mu X]\) uses Equation 19.3 — \(\mu\) is a constant, not random.

4. Two variables have \(\sigma_X = 2\), \(\sigma_Y = 3\), \(\rho = 0.5\). Find \(\operatorname{Cov}(X,Y)\) and \(\operatorname{Var}(X+Y)\).

Solution

\(\operatorname{Cov} = \rho\sigma_X\sigma_Y = 0.5 \times 2 \times 3 = 3\), so \(\operatorname{Var}(X+Y) = 4 + 9 + 2(3) = 19\).

cv <- 0.5 * 2 * 3
c(covariance = cv, var_sum = 2^2 + 3^2 + 2 * cv)

covariance    var_sum 
         3         19

Had they been uncorrelated the answer would be 13. Positive correlation makes a sum more variable, which is the whole argument for diversification: uncorrelated or negatively correlated components give a portfolio less variance than the sum of its parts.

5. Verify that \(\mathbb{E}[X+Y] = \mathbb{E}[X] + \mathbb{E}[Y]\) even for strongly dependent variables.

Solution

xdep <- rnorm(5000)
ydep <- xdep^2 # completely determined by xdep
c(
  E_sum = mean(xdep + ydep),
  sum_of_E = mean(xdep) + mean(ydep)
)

    E_sum  sum_of_E 
0.9815641 0.9815641

Identical, as they must be. Linearity of expectation never requires independence — unlike variance, where the covariance term appears. This is why expected-value arguments work in situations where variance calculations become intractable.

6. For the joint table above, compute \(\mathbb{E}[Y \mid X = 1]\) and verify the law of total expectation.

Solution

cond_y_given_x <- joint / rowSums(joint)
round(cond_y_given_x, 4)

   Y
X        0      1
  0 0.3333 0.6667
  1 0.4286 0.5714

Ey_given <- as.vector(cond_y_given_x %*% lvl)
c(
  E_Y_given_X0 = Ey_given[1],
  E_Y_given_X1 = Ey_given[2],
  tower = sum(Ey_given * rowSums(joint)),
  E_Y = EYj
)

E_Y_given_X0 E_Y_given_X1        tower          E_Y 
   0.6666667    0.5714286    0.6000000    0.6000000

Each row of the conditional table sums to 1 — that is the renormalization in Equation 18.2. Averaging the two conditional means, weighted by \(P(X=x)\), recovers \(\mathbb{E}[Y]\) exactly, which is Equation 19.9.