4  Random Variables

Previously, we computed a distribution given the data, whereas now we generate individual data points given the distribution. This chapter introduces random variables and the two basic kinds of sample space (discrete and continuous), along with the most common distributions in each family.

To talk about data we have not yet observed, we need an object that names the set of possible outcomes and the long-run frequency of each.

ImportantKey Definition

A random variable \(X_{i}\) is a variable whose value is generated according to a probability distribution. Its sample space is the set of all values it can take, and its probability assigns a long-run frequency to each subset of the sample space. A realized value \(\hat{x}_{i}\) is a specific observed value drawn from the random variable’s distribution.

The random-variable framing is useful for reasoning about a data generating process before any data are observed. For example, consider flipping a coin before knowing whether it lands on heads or tails: we denote the unflipped coin as the random variable \(X_{i}\), with sample space \(\{0, 1\}\) (\(x=1\) is heads, \(x=0\) is tails) and probability \(1/2\) on each outcome. A specific flipped coin that has landed heads or tails is denoted \(\hat{x}_{i}\), mirroring the hat notation from the previous chapter.

There are two basic types of sample spaces: discrete (encompassing cardinal-discrete, factor-ordered, and factor-unordered data) and continuous. This leads to two types of random variables: discrete and continuous. However, each type has many different probability distributions. The most common random variables easily accessible in R.

Cumulative Probability

The first job of any random-variable description is to assign a probability to events of the form “the value falls at or below \(x\)”.

ImportantKey Definition

The cumulative distribution function (CDF) of a random variable \(X_{i}\) is \[F(x) = \text{Prob}(X_{i} \leq x),\] the theoretical probability that \(X_{i}\) takes a value at or below \(x\). The CDF is the theoretical analogue of the empirical CDF \(\hat{F}\) from the previous chapter.

You can think of \(F(x)\) as the ECDF for a dataset with an infinite number of observations; equivalently, the ECDF is an empirical version of the CDF applied to observed data. The CDF is useful because it packages every probability of interest into a single function, and works for both continuous random variables and ordered factors.

Probability Rules

After introducing different random variables, we will also cover some basic implications of their CDF. Intuitively, probabilities must sum up to one. So we can compute \(\text{Prob}(X_{i} > x) = 1- F(x)\). We also have two “in” and “out” probabilities.

The probability of \(X_{i}\leq b\) and \(X_{i}> a\) can be written in terms of falling into a range \[\text{Prob}(X_{i}~\text{in}~(a,b])=\text{Prob}(a < X_{i} \leq b) = F(b) - F(a).\] The opposite probability is \[\begin{aligned} \text{Prob}(X_{i} \text{ out of } (a,b]) &= \text{Prob}(X_{i} \leq a \text{ or } X_{i} > b) = F(a) + [1- F(b)] \\ &= 1 - [F(b) - F(a)] = 1 - \text{Prob}( X_{i}~\text{in}~(a,b]). \end{aligned}\]

The endpoints matter when \(X_{i}\) is discrete, because a single value can have positive probability. Including the left endpoint adds the probability of exactly \(a\), so \(\text{Prob}(X_{i}~\text{in}~[a,b]) = F(b) - F(a) + \text{Prob}(X_{i}=a)\). The endpoints do not matter when \(X_{i}\) is continuous (no single value has positive probability).

4.1 Discrete Random Variables

A discrete random variable can take one of several values in a set. E.g., any number in \(\{1,2,3,...\}\). With a discrete random variable, individual values can have positive probability, so we can describe the distribution by listing the probability of each value.

ImportantKey Definition

The probability mass function (PMF) of a discrete random variable \(X_{i}\) is \[p(x) = \text{Prob}(X_{i}=x),\] the theoretical probability that \(X_{i}\) takes exactly the value \(x\). The PMF is the theoretical analogue of the empirical mass function \(\hat{p}(x)\).

You can think of \(p(x)\) as a proportions bar plot for a dataset with an infinite number of observations; equivalently, the bar plot is an empirical version of the PMF applied to observed data. The CDF accumulates the mass function, \(F(x)=\sum_{x' \leq x} p(x')\), and conversely the probability of any value is the size of the jump in \(F\) at that value. The PMF is useful for reading off the probability of a particular outcome, such as rolling exactly a \(3\), and it is how you specify a discrete random variable in R.

In R, you can simulate random discrete data using the sample function. To do so, you must specify each outcome \(x\) and its probability \(p(x)\).

Code 
x <- c(0, 1)
x_probs <- c(1/2, 1/2)
sample(x, 1, prob=x_probs, replace=TRUE) # One Flip
## [1] 0
sample(x, 1, prob=x_probs, replace=TRUE) # Another Flip
## [1] 1
sample(x, 4, prob=x_probs, replace=TRUE) # 4 Flips 
## [1] 0 0 0 0

Bernoulli

The simplest random variable models a single yes/no event, like one flip of a coin. A Bernoulli random variable takes two values, \(x=1\) (success/heads) with probability \(\theta\) and \(x=0\) (failure/tails) with probability \(1-\theta\); the single parameter \(\theta~\text{in}~[0,1]\) controls the probability of success.

Bernoulli variables are useful whenever an outcome is binary (heads or tails, success or failure, present or absent), and are also the building block from which richer discrete distributions are assembled out of independent yes/no events. Writing the distribution out in full helps fix the notation we will reuse for richer cases: \[\begin{aligned} X_{i} &~\text{in}~\{0,1\} \\ p(x) &= \begin{cases} 1-\theta & x=0 \\ \theta & x=1 \\ 0 & \text{Otherwise} \end{cases} \\ F(x) &= \begin{cases} 0 & x<0 \\ 1-\theta & x~\text{in}~[0,1) \\ 1 & x\geq 1 \end{cases} \end{aligned}\]

The theoretical mean of a Bernoulli variable works out to \(\theta\).

Here is an example of the Bernoulli distribution. While you might get all heads (or all tails) in the first few coin flips, the ratios level out to their theoretical values after many flips.

Code 
x <- c(0, 1)
x_probs <- c(3/4, 1/4)

sample(x, 1, prob=x_probs, replace=TRUE) # 1 Flip
## [1] 0
x_sim <- sample(x, 400, prob=x_probs, replace=TRUE) #400 Flips

# Plot proportions
proportions <- table(x_sim)/length(x_sim)
plot(proportions, 
    col=grey(0, 0.5),
    xlab='Flip Outcome', 
    ylab='Proportion',
    main=NA,
    ylim=c(0,0.8))
points(c(0, 1), c(.75, .25), pch=16, col=rgb(0, 0, 1, .8)) # Theoretical values

Code 

# Plot CDF
plot( ecdf(x_sim), col=grey(0, .5),
    pch=16, main=NA) #Empirical

Code 
#points(c(0, 1), c(.75, 1), pch=16, col='blue') # Theoretical

Discrete Uniform

Discrete numbers with equal probability, such as a die with \(K\) sides. \[\begin{aligned} X_{i} &~\text{in}~\{1,...K\} \\ p(x) &= \begin{cases} 1/K & x~\text{in}~\{1,...K\} \\ 0 & \text{Otherwise} \end{cases} \\ F(x) &= \begin{cases} 0 & x<1 \\ 1/K & x~\text{in}~[1,2) \\ 2/K & x~\text{in}~[2,3) \\ \vdots & \\ 1 & x\geq K \end{cases} \end{aligned}\]

Here is an example with \(K=4\). E.g., rolling a four-sided die.

The probability of a value smaller than or equal to \(3\) is \(\text{Prob}(X_{i} \leq 3)=1/4 + 1/4 + 1/4 = 3/4\).

The probability of a value larger than \(3\) is \(\text{Prob}(X_{i} > 3) = 1-\text{Prob}(X_{i} \leq 3)=1/4\).

The probability of a value \(>\) 1 and \(\leq 3\) is \(\text{Prob}(1 < X_{i} \leq 3) = \text{Prob}(X_{i} \leq 3) - \text{Prob}(X_{i} \leq 1) = 3/4 - 1/4 = 2/4\).1

Including \(1\) changes the answer, since the event now also contains rolling a \(1\): \(\text{Prob}(1 \leq X_{i} \leq 3) = \text{Prob}(X_{i} \leq 3) - \text{Prob}(X_{i} \leq 1) + \text{Prob}(X_{i}=1) = 3/4 - 1/4 + 1/4 = 3/4\).

The probability of a value \(\leq\) 1 or \(> 3\) is \(\text{Prob}(X_{i} \leq 1 \text{ or } X_{i} > 3) = \text{Prob}(X_{i} \leq 1) + \left[ 1- \text{Prob}(X_{i} \leq 3) \right] = 1/4 + [1 - 3/4]=2/4\).

Code 
x <- c(1, 2, 3, 4)
x_probs <- c(1/4, 1/4, 1/4, 1/4)
# sample(x, 1, prob=x_probs, replace=T) # 1 roll
x1_hat <- sample(x, 2000, prob=x_probs, replace=TRUE) # 2000 rolls

# Plot Long run proportions
proportions <- table(x1_hat)/length(x1_hat)
plot(proportions, col=grey(0, .5),
    xlab='Outcome', ylab='Proportion', main=NA)
points(x, x_probs, pch=16, col=rgb(0, 0, 1, .8)) # Theoretical values

Code 

# Hist w/ Theoretical Counts
# hist(x1_hat, breaks=50, border=NA, main=NA, ylab='Count')
# points(x, x_probs*length(x1_hat), pch='-')

# Alternative Plot
plot( ecdf(x1_hat), pch=16, col=grey(0, .5), main=NA)

Code 

# Alternative Plot 2
#props <- table(x1_hat)
#barplot(props, ylim = c(0, 0.35), ylab = 'Proportion', xlab = 'Value')
#abline(h = 1/4, lty = 2)

Note that the Discrete Uniform distribution generalizes to arbitrary sets, such as \(\{0.0, 0.25, 0.5, 0.75, 1.0\}\) instead of \(\{1,2,3,4,5\}\). We will not exploit the generalization in this class.

Multinoulli

Numbers \(1,...K\), with each outcome \(k\) having probability \(\theta_{k}\) and satisfying \(\theta_{1} + \theta_{2} + ... = 1\). \[\begin{aligned} X_{i} &~\text{in}~\{1,...K\} \\ p(x) &= \begin{cases} \theta_{1} & x=1 \\ \theta_{2} & x=2 \\ \vdots & \\ \theta_{K} & x=K \\ 0 & \text{Otherwise} \end{cases} \\ F(x) &= \begin{cases} 0 & x<1 \\ \theta_{1} & x~\text{in}~[1,2) \\ \theta_{1} + \theta_{2} & x~\text{in}~[2,3) \\ \vdots & \\ 1 & x\geq K \end{cases} \end{aligned}\]

Here is an empirical example with three outcomes

Code 
x <- c(1, 2, 3)
x_probs <- c(3/10, 1/10, 6/10)
sum(x_probs)
## [1] 1
x2_hat <- sample(x, 2000, prob=x_probs, replace=TRUE) # sample of 2000

# Plot Long run proportions
proportions <- table(x2_hat)/length(x2_hat)
plot(proportions, col=grey(0, .5), ylim=c(0, 0.7),
    xlab='Outcome', ylab='Proportion', main=NA)
points(seq(x), x_probs, pch=16, col=rgb(0, 0, 1, .8)) # Theoretical values

Code 

# Histogram version
# x2_alt_hat <- x1_hat
# x2_alt_hat[x2_alt_hat=='A'] <- 1
#x2_alt_hat[x2_alt_hat=='B'] <- 2
#x2_alt_hat[x2_alt_hat=='C'] <- 3
#x2_alt_hat <- as.numeric(x1_alt_hat)
#hist(x2_alt_hat, breaks=50, border=NA,
#    main=NA, ylab='Count')
#points(x, x_probs*length(x2_alt_hat), pch=16) ## Theoretical Counts

# Alternative Plot
# plot( ecdf(x2_hat), pch=16, col=grey(0, .5), main=NA)

We can also replace numbers with letters \((A,...Z)\) or names \((\text{John}, \text{Jamie}, ...)\) although we must be careful with the CDF when there is no longer a natural ordering. When the outcomes are not cardinal data, the Multinoulli random variable is typically referred to as a Categorical random variable. The CDF \(F(x)=\text{Prob}(X_{i} \leq x)\) requires an ordering on the outcomes, and with labels like “my car dies” and “it rains next Tuesday”, you do not have any inherent order and the CDF is an arbitrary artifact.

Suppose there are three possible outcomes of an experiment, \(\{\text{my car dies}, \text{it rains next Tuesday},\text{a cat is born}\}\), which have corresponding probabilities \(\{3/10, 1/10, 6/10 \}\) that are known theoretically. Compute the probability that my car dies or a cat is born.

Draw Random Letters with increasing probability

Code 
x <- LETTERS
x_probs <- seq(1, length(x))
x_probs <- x_probs/sum(x_probs) # probs must sum to 1
x_sim <- sample(x, 1, prob=x_probs, replace=TRUE)
x_sim
## [1] "M"

# Draw Random Letters 1000 Times
x_sim <- sample(x, 1000, prob=x_probs, replace=TRUE)
x_sim <- factor(unlist(x_sim), levels=LETTERS, ordered=FALSE)

4.2 Continuous Random Variables

A continuous random variable can take one value out of an uncountably infinite number. E.g., any number between \(0\) and \(1\) with any number of decimal points. With a continuous random variable, the probability of any individual point is zero, so we cannot describe the distribution by listing \(\text{Prob}(X_{i} = x)\) for each \(x\); instead, we describe it with the CDF or with a density.

ImportantKey Definition

The probability density function (PDF) \(f(x)\) of a continuous random variable describes how concentrated the probability is near each value of \(x\). The probability that \(X_{i}\) falls in an interval \([a,b]\) is the area under \(f\) from \(a\) to \(b\), equivalently \(F(b) - F(a)\). \(f(x)\) itself is not a probability; only its area over an interval is.

The PDF is useful for intuition: probabilities are pictured as areas, peaks of \(f\) mark where outcomes are likely, and the shape of \(f\) tells you whether the distribution is symmetric, skewed, or heavy-tailed. The CDF is the workhorse for numerical calculation, since \(F(x)\) already encodes the cumulative area under \(f\) from \(-\infty\) to \(x\): where \(F\) rises steeply \(f\) is large, and where \(F\) is flat \(f\) is near zero. Just as \(F\) can be thought of as the ECDF \(\hat{F}\) with an infinite amount of data, \(f\) can be thought of as a histogram \(\hat{f}\) with an infinite amount of data; equivalently, the histogram is an empirical version of the PDF applied to observed data.

In R, you work with continuous distributions using these functions

  • dXXX density, \(f(x)\) for particular values \(x\)
  • pXXX distribution, \(F(x)\) for particular values \(x\)
  • rXXX generate random variables according to the distribution

where XXX is the distribution name (unif, exp, norm)

Continuous Uniform

The simplest continuous random variable is the one with no preferred region: every value in an interval is equally likely. A Continuous Uniform random variable on the interval \([\theta_{L},\theta_{H}]\) is equally likely to land anywhere in \([\theta_{L},\theta_{H}]\), with density \(f(x) = 1/(\theta_{H}-\theta_{L})\) for \(x~\text{in}~[\theta_{L},\theta_{H}]\) and \(0\) otherwise.

The Continuous Uniform is useful as a reference distribution: a flat density makes it a natural “no information” baseline, and R’s runif is the foundation underlying many other random-number routines. The full specification gathers sample space, density, and cumulative function in one place: \[\begin{aligned} X_{i} &~\text{in}~[\theta_{L},\theta_{H}] \\ f(x) &= \begin{cases} 1/(\theta_{H}-\theta_{L}) & x~\text{in}~[\theta_{L},\theta_{H}] \\ 0 & \text{Otherwise} \end{cases} \\ F(x) &= \begin{cases} 0 & x < \theta_{L} \\ (x-\theta_{L})/(\theta_{H}-\theta_{L}) & x~\text{in}~[\theta_{L},\theta_{H}] \\ 1 & x > \theta_{H}. \end{cases} \end{aligned}\]

The density \(f\) and the CDF \(F\) hold the same information, and you can move between them. Consider a Continuous Uniform on \([-1, 2]\), so its density is flat: \(f(x)=1/3\) for \(x~\text{in}~[-1, 2]\).

To build \(F\) from \(f\), accumulate area. The area under \(f\) from \(-1\) to \(x\) is a rectangle of width \((x-(-1))=x+1\) and height \(1/3\), so \(F(x)=(x+1)/3\). Check: \(F(-1)=0\), \(F(2)=1\), and \(F(0.5)=1.5/3=0.5\).

To recover \(f\) from \(F\), take the slope. The CDF \(F(x)=(x+1)/3\) is a straight line of slope \(1/3\), so \(f(x)=1/3\).

Code 
x <- seq(-1, 2, by=0.01)
par(mfrow=c(1, 2))
# Density f(x): flat at 1/3
plot(x, dunif(x, -1, 2), type='l', col=rgb(0, 0, 1, .8),
    ylim=c(0, 0.5), xlab='x', ylab='f(x)', main=NA)
title('Density f(x)', font.main=1)
# CDF F(x): accumulated area under f
plot(x, punif(x, -1, 2), type='l', col=rgb(0, 0, 1, .8),
    xlab='x', ylab='F(x)', main=NA)
title('CDF F(x)', font.main=1)

Note that the standard uniform uses the unit interval \([0,1]\) and is the special case \(\theta_{L}=0, \theta_{H}=1\). The density is \(f(x)=1\) and CDF \(F(x)=x\).

For the unit-interval special case (\(\theta_{L}=0, \theta_{H}=1\)), the probability of a value being exactly \(0.25\) is \(\text{Prob}(X_{i} =0.25)=0\).

Since \(\text{Prob}(X_{i}=a)=0\), the endpoints do not matter for a continuous random variable: \([a,b]\), \((a,b]\), \([a,b)\), and \((a,b)\) all have probability \(F(b)-F(a)\). This differs from the four-sided die, where \(\text{Prob}(1 \leq X_{i} \leq 3) = 3/4\) but \(\text{Prob}(1 < X_{i} \leq 3) = 2/4\).

The probability of a value smaller than \(0.25\) is \(F(0.25)=0.25\).

The probability of a value larger than \(0.25\) is \(1-F(0.25)=0.75\).

The probability of a value in \((0.25,0.75]\) is \(\text{Prob}(0.25 < X_{i} \leq 0.75) = \text{Prob}(X_{i} \leq 0.75) - \text{Prob}(X_{i} \leq 0.25) = 0.75 - 0.25 = 0.5\).

The probability of a value in \((0.2,0.7]\) is \(\text{Prob}(0.2 < X_{i} \leq 0.7) = \text{Prob}(X_{i} \leq 0.7) - \text{Prob}(X_{i} \leq 0.2) = 0.7 - 0.2 = 0.5\).

The probability of a value outside of \((0.2,0.7]\) is \(\text{Prob}(X_{i} \leq 0.2 \text{ or } x > 0.7) = 0.2 + [1-0.7]=0.5\). Alternatively, you can compute \(1- \text{Prob}(0.2 < X_{i} \leq 0.7)=1-0.5=0.5\).

Code 
# Prob. < 0.25
punif(0.25)
## [1] 0.25

# Prob. > 0.25
1-punif(0.25)
## [1] 0.75

# Prob. in (0.2, 0.7]
F_two <- punif( c(0.2, 0.7) )
F_in <- F_two[2] - F_two[1]
F_in
## [1] 0.5

# Prob. out of (0.2, 0.7]
F_out <- 1- F_in
F_out
## [1] 0.5
Code 
runif(3) # 3 draws
## [1] 0.70842381 0.03178858 0.02837421

# Empirical Density 
x3_hat <- runif(2000)
hist(x3_hat, breaks=20, border=NA, main=NA, freq=FALSE)
# Theoretical Density
x <- seq(-0.1, 1.1, by=.001)
fx <- dunif(x)
lines(x, fx, col=rgb(0, 0, 1, .8))

Code 

# CDF example 1
P_low <- punif(0.25)
P_low
## [1] 0.25
# Uncomment to show via PDF
# x_low <- seq(0, 0.25, by=.001)
# fx_low <- dunif(x_low)
# polygon( c(x_low, rev(x_low)), c(fx_low, fx_low*0),
#    col=rgb(0, 0, 1, .25), border=NA)
    
# CDF example 2
P_high <- 1-punif(0.25)
P_high
## [1] 0.75
# Uncomment to show via PDF
# x_high <- seq(0.25, 1, by=.001)
# fx_high <- dunif(x_high)
# polygon( c(x_high, rev(x_high)), c(fx_high, fx_high*0),
#    col=rgb(0, 0, 1, .25), border=NA)
    
# CDF example 3
P_mid <- punif(0.75) - punif(0.25)
P_mid
## [1] 0.5
# Uncomment to show via PDF
# x_mid <-  seq(0.25, 0.75, by=.001)
# fx_mid <- dunif(x_mid)
# polygon( c(x_mid, rev(x_mid)), c(fx_mid, fx_mid*0),
#    col=rgb(0, 0, 1, .25), border=NA)

Suppose \(X_{i}\) is a random variable continuously distributed over \(\theta_{L}=-2\) and \(\theta_{H}=2\). What is the probability of a value smaller than \(0.25\)? First use the computer to suggest an answer: simulate \(1000\) draws and then make a histogram and an ECDF. Then find the answer mathematically using the CDF. Finally, verify the answer is intuitively correct in a figure of the PDF. You should draw by hand both the CDF and the PDF with correct axes labels and marking clearly the probability of a value smaller than \(0.25\).

Code 
# Simulation
x_sim <- runif(1000, -2, 2)
# ... 

# Draw by computer the CDF, F(x)
x <- seq(-3, 3, by=0.005)
Fx <- punif(x, -2, 2)
plot(x, Fx, type='l', col=rgb(0, 0, 1, .8))
# Answer
Fstar <- punif(0.25, -2, 2)
# Visualize Answer
segments(0.25, 0, 0.25, Fstar,
    col=rgb(0, 0, 1, 0.25))

Code 

# Draw by computer the PDF, f(x)
fx <- dunif(x, -2, 2)
plot(x, fx, type='l', col=rgb(0, 0, 1, .8))
# Visualize Answer
fstar <- dunif(0.25, -2, 2)
polygon(
    c(-2, 0.25, 0.25, -2), 
    c(0, 0, fstar, fstar),
    col=rgb(0, 0, 1, 0.25), border=NA)

Suppose the flight time between Calgary and Kamloops is Uniformly distributed between \(68\) and \(78\) minutes. According to Air Canada the flight takes \(70\) minutes. What is the probability that the flight will be late?

Exponential

The sample space is any positive number.2 An Exponential random variable has a single rate parameter, \(\lambda>0\). \[\begin{aligned} X_{i} &~\text{in}~[0,\infty) \\ f(x) &= \lambda \exp\left\{ -\lambda x \right\} \\ F(x) &= \begin{cases} 0 & x < 0 \\ 1- \exp\left\{ -\lambda x \right\} & x \geq 0. \end{cases} \end{aligned}\]

Code 
rexp(3) # 3 draws
## [1] 3.437327 1.634091 0.608449

x5_hat <- rexp(2000)
hist(x5_hat, breaks=20,
    border=NA, main=NA,
    freq=FALSE, ylim=c(0, 1), xlim=c(0, 10))

x <- seq(0, 10, by=.1)
fx <- dexp(x)
lines(x, fx, col=rgb(0, 0, 1, .8))

Suppose the lifetime of a battery is an exponential random variable with \(\lambda=1/50\). Using the computer, find the probability that the lifetime is \(< 10\) hours. Find the probability that the lifetime is \(\geq 100\) hours. Use the computer to find the probability that the lifetime is between \(10\) and \(100\) hours.

Code 
pexp(10, 1/50)
## [1] 0.1812692

Note that there is a discrete version of the exponential, called the geometric distribution.

Code 
y_hat <- rexp(5000, rate=.2)

x_sim <- rgeom(5000, 1-exp(-.2))

par(mfrow=c(1, 2))
bks <- seq(0, 60, by=0.25)   # span both samples
hist(x_sim, breaks=bks, freq=FALSE, border=NA, main=NA, xlab='x')
title('Geometric', font.main=1)
hist(floor(y_hat), breaks=bks, freq=FALSE, border=NA, main=NA, xlab='x')
title('Exponential (floored)', font.main=1)

Normal (Gaussian)

For continuous data spread symmetrically around a center, the default model is the bell curve. A Normal (or Gaussian) random variable has a bell-shaped density on the real line, with two parameters: the mean \(\mu\) (the center) and the standard deviation \(\sigma\) (the spread). The Normal is useful as the default model for measurements that cluster around a typical value with symmetric scatter (heights, test scores, measurement errors), and is sometimes called the Gaussian distribution. We call it “Normal” because it appears again and again, both as a description of real data and as the limiting shape for sample averages. The full PDF spells out the parameters \(\mu~\text{in}~(-\infty,\infty)\) and \(\sigma > 0\): \[\begin{aligned} X_{i} &~\text{in}~(-\infty,\infty) \\ f(x) &= \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left\{ \frac{-(x-\mu)^2}{2\sigma^2} \right\} \end{aligned}\] The CDF, \(F(x)\), has no closed expression. Despite the complicated formula, the shape is simple: a bell centered at \(\mu\), with most of the probability within a few multiples of \(\sigma\) on either side.

Code 
rnorm(3) # 3 draws
## [1] -1.7749409  0.3200881 -1.1568385

x6_hat <- rnorm(2000)
hist(x6_hat, breaks=20,
    border=NA, main=NA,
    freq=FALSE, ylim=c(0, .4), xlim=c(-4, 4))

x <- seq(-10, 10, by=.025)
fx <- dnorm(x)
lines(x, fx, col=rgb(0, 0, 1, .8))

Even thought the distribution function is complex, we can compute CDF values using the computer.

Code 
pnorm( c(-1.645, 1.645) ) # 10%
## [1] 0.04998491 0.95001509
pnorm( c(-2.576, 2.576) ) #  1%
## [1] 0.004997532 0.995002468

Suppose \(X_{i}\) is a random variable with a normal distribution with \(\mu=0\) and \(\sigma=1\). Intuitively depict \(\text{Prob}(X_{i}~\text{in}~[0.2, 0.8])\) by drawing an area under the density function. Numerically estimate that same probability using the CDF.

Code 
plot( ecdf(x6_hat), main=NA) # Empirical

x <- seq(-4, 4, by=.01) # Theoretical
Fx <- pnorm(x, 0, 1)
lines(x, Fx, col=rgb(0, 0, 1, .8))

# Middle Interval Example
F2 <- pnorm(0.2, 0, 1)
F8 <- pnorm(0.8, 0, 1)
F_2_8 <- F8 - F2
F_2_8
## [1] 0.2088849

# Visualize
title('Middle between 0.2 and 0.8', font.main=1)
segments( 0.2, F2, -5, F2, col=rgb(1, 0, 0, .8))
segments( 0.8, F8, -5, F8, col=rgb(1, 0, 0, .8))

Suppose that your health status is a normally distributed random variable with \(\mu=2\) and \(\sigma=3\). If we randomly sample one person, what is the probability there health status is higher than \(4\)?

Code 
# Start with a simulation of 1000 people to build intuition
x_sim <- rnorm(1000, 2, 3)
hist(x_sim, freq=FALSE, border=NA, main=NA)

Code 
sum(x_sim>4)/1000
## [1] 0.257

# Do an exact calculation
1-pnorm(4, 2, 3)
## [1] 0.2524925

Draw the Middle \(90\%\) of a normal distribution

Code 
# PDF
x <- seq(-10, 10, by=.025)
fx <- dnorm(x)
plot(x, fx,
    col=rgb(0, 0, 1, .8), type='l',
    ylab='Theoretical Density: f(x)',
    main=NA)
title('Middle 90%', font.main=1)
# Show Middle 90%
x_90 <- seq(-1.645, 1.645, by=.025)
fx_90 <- dnorm(x_90)
polygon( c(x_90, rev(x_90)), c(fx_90, fx_90*0),
    col=rgb(0, 0, 1, .25), border=NA)

Code 
pnorm(1.645)-pnorm(-1.645)
## [1] 0.9000302

Suppose scores in a math class are approximately normally distributed with \(\mu=70\) and \(\sigma=10\). If you select one student at random, what is the probability their score is higher than \(90\)? Is \(\text{Prob}(X_{i}\geq 90)\) higher if \(\mu=80, \sigma=10\)? What about \(\mu=70, \sigma=20\)?

Code 
# Case 1: baseline, mu=70 and sigma=10
1 - pnorm(90, mean=70, sd=10)
## [1] 0.02275013

# Case 2: shift the mean up to 80, same sigma
1 - pnorm(90, mean=80, sd=10)
## [1] 0.1586553

# Case 3: keep mean at 70, double the spread to sigma=20
1 - pnorm(90, mean=70, sd=20)
## [1] 0.1586553

Raising \(\mu\) from \(70\) to \(80\) moves the bell closer to \(90\), so the right-tail probability rises sharply. Doubling \(\sigma\) at fixed \(\mu\) also raises the right-tail probability because the bell flattens and more mass spreads into the tails.

4.3 Exercises

  1. Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.

  2. A Multinoulli random variable has outcomes \(\{1,2,3,4\}\) with probabilities \(\{0.1, 0.3, 0.4, 0.2\}\). Using the CDF, compute \(\text{Prob}(X_{i} \leq 2)\), \(\text{Prob}(X_{i} > 3)\), and \(\text{Prob}(1 < X_{i} \leq 3)\). Show your work.

  3. Suppose \(X_{i}\) is normally distributed with \(\mu=10\) and \(\sigma=4\). Using pnorm, compute the probability that \(X_{i}\) falls between \(6\) and \(14\). Then compute the probability that \(X_{i}\) is more than two standard deviations away from the mean, i.e., \(\text{Prob}(X_{i} < 2 \text{ or } X_{i} > 18)\).

  4. In R, generate \(2000\) draws from an exponential distribution with \(\lambda=0.5\) using rexp(2000, rate=0.5). Plot a histogram with freq=FALSE and overlay the theoretical density using lines and dexp. Then use pexp to compute the exact probability that a draw exceeds \(4\).

Further Reading

Many random variables are related to each other. Also note that numbers randomly generated on your computer cannot be truly random, they are “pseudorandom.”

Recall

This chapter introduced random variables and distinguished discrete families (Bernoulli, Uniform, Multinoulli) from continuous families (Uniform, Exponential, Normal), describing each through the CDF \(F(x)=\text{Prob}(X_{i} \leq x)\) and (for continuous variables) the PDF \(f(x)\). The four-sided die example tied the CDF rules together: with \(\text{Prob}(X_{i} \leq 3) = 3/4\), the right-tail probability is \(\text{Prob}(X_{i} > 3) = 1/4\), the in-interval probability is \(\text{Prob}(1 < X_{i} \leq 3) = 3/4 - 1/4 = 2/4\), and the out-of-interval probability is \(\text{Prob}(X_{i} \leq 1 \text{ or } X_{i} > 3) = 1/4 + 1/4 = 2/4\).

  1. This is the general formula using CDFs, and you can verify it works in this instance by directly adding the probability of each 2 or 3 event: \(\text{Prob}(X_{i} = 2) + \text{Prob}(X_{i} = 3) = 1/4 + 1/4 = 2/4\).↩︎

  2. In other classes, you may further distinguish types of random variables based on whether their maximum value is theoretically finite or infinite.↩︎