SD04 — Random Variables
Course: Statistics and Probability for Data Science (DS14270).
Source deck: Samer El Kababji, SD04 — Random Variables.
This lesson develops the deck’s definitions of expectation, variance, standard deviation, percentiles, and quantiles for discrete and continuous random variables. It reproduces the machine-breakdown and exponential worked examples in executable R.
The expected value is a probability-weighted average of the possible values of a random variable.
For a discrete random variable \(X\) with possible values \(D\) and PMF \(p\),
\[ E(X)=\mu_X=\sum_{x\in D}x p(x). \]
discrete_expectation <- function(x, probability) {
stopifnot(length(x) == length(probability),
all(probability >= 0),
isTRUE(all.equal(sum(probability), 1)))
sum(x * probability)
}
discrete_expectation(
x = c(0, 1, 2),
probability = c(0.5, 1 / 6, 1 / 3)
)[1] 0.8333333For a continuous random variable \(X\) with PDF \(f\),
\[ E(X)=\mu_X=\int_{-\infty}^{\infty}x f(x)\,dx. \]
continuous_expectation <- function(density, lower = -Inf, upper = Inf) {
integrate(function(x) x * density(x), lower, upper)$value
}
unit_interval_density <- function(x) {
ifelse(x >= 0 & x <= 1, 1, 0)
}
continuous_expectation(unit_interval_density, 0, 1)[1] 0.5The expectation describes the long-run probability center. It need not be a possible value: the expected result of a fair six-sided die is \(3.5\).
Variance is the expected squared distance from the mean. For a discrete random variable,
\[ \operatorname{Var}(X) =\sum_{x\in D}(x-\mu)^2p(x) =E(X^2)-\{E(X)\}^2. \]
For a continuous random variable,
\[ \operatorname{Var}(X) =\int_{-\infty}^{\infty}(x-\mu)^2f(x)\,dx =E(X^2)-\{E(X)\}^2. \]
The standard deviation returns to the units of \(X\):
\[ \sigma_X=\sqrt{\operatorname{Var}(X)}. \]
discrete_moments <- function(x, probability) {
mean_x <- discrete_expectation(x, probability)
second_moment <- sum(x^2 * probability)
variance_x <- second_moment - mean_x^2
c(
mean = mean_x,
second_moment = second_moment,
variance = variance_x,
standard_deviation = sqrt(variance_x)
)
}
discrete_moments(1:6, rep(1 / 6, 6)) mean second_moment variance standard_deviation
3.500000 15.166667 2.916667 1.707825 For \(0<p<1\), the \(p\) quantile is defined by the generalized inverse of the CDF:
\[ Q_p=\inf\{x:F(x)\geq p\}. \]
This definition handles both smooth continuous CDFs and jumping discrete CDFs. The median is \(Q_{0.50}\); the three commonly reported quartiles are
\[ Q_1=Q_{0.25}, \qquad Q_2=Q_{0.50}, \qquad Q_3=Q_{0.75}. \]
The interquartile range is
\[ \operatorname{IQR}=Q_3-Q_1. \]
For a discrete PMF, the first support value whose cumulative probability reaches \(p\) implements the generalized inverse.
discrete_quantile <- function(x, probability, p) {
stopifnot(length(x) == length(probability),
all(diff(x) > 0),
isTRUE(all.equal(sum(probability), 1)),
all(p > 0 & p < 1))
cumulative_probability <- cumsum(probability)
vapply(p, function(level) {
x[which(cumulative_probability >= level)[1]]
}, numeric(1))
}
discrete_quantile(
x = c(1, 2, 4, 8, 16),
probability = c(0.05, 0.10, 0.35, 0.40, 0.10),
p = c(0.25, 0.50, 0.75)
)[1] 4 4 8Let \(X\) be a machine-breakdown cost with PMF
| Cost \(x\) | 50 | 200 | 350 |
|---|---|---|---|
| \(P(X=x)\) | 0.3 | 0.2 | 0.5 |
First,
\[ P(X\geq200)=P(X=200)+P(X=350)=0.2+0.5=0.7. \]
breakdown_cost <- c(50, 200, 350)
breakdown_probability <- c(0.3, 0.2, 0.5)
P_at_least_200 <- sum(
breakdown_probability[breakdown_cost >= 200]
)
P_at_least_200[1] 0.7The mean is
\[ E(X)=50(0.3)+200(0.2)+350(0.5)=230. \]
The second moment and variance are
\[ E(X^2)=50^2(0.3)+200^2(0.2)+350^2(0.5)=70{,}000, \]
\[ \operatorname{Var}(X)=70{,}000-230^2=17{,}100, \]
so
\[ \operatorname{SD}(X)=\sqrt{17{,}100}\approx130.77. \]
breakdown_moments <- discrete_moments(
breakdown_cost,
breakdown_probability
)
breakdown_moments mean second_moment variance standard_deviation
230.000 70000.000 17100.000 130.767 stopifnot(isTRUE(all.equal(
unname(breakdown_moments[c("mean", "second_moment", "variance")]),
c(230, 70000, 17100)
)))The CDF values on the support are
\[ F(50)=0.3, \qquad F(200)=0.5, \qquad F(350)=1. \]
Using \(Q_p=\inf\{x:F(x)\geq p\}\),
\[ Q_1=50, \qquad Q_2=200, \qquad Q_3=350, \qquad Q_{0.40}=200. \]
probability_levels <- c(Q1 = 0.25, Q2 = 0.50, Q3 = 0.75,
Q40 = 0.40)
breakdown_quantiles <- discrete_quantile(
breakdown_cost,
breakdown_probability,
probability_levels
)
c(
breakdown_quantiles,
IQR = breakdown_quantiles["Q3"] - breakdown_quantiles["Q1"]
) Q1 Q2 Q3 Q40 IQR.Q3
50 200 350 200 300 These are distribution quantiles, obtained from the PMF. They should not be confused with an interpolated sample-quantile convention applied to a data vector.
Let \(X\) have an exponential distribution with rate \(\lambda>0\):
\[ f(x)= \begin{cases} \lambda e^{-\lambda x}, & x\geq0,\\ 0, & x<0. \end{cases} \]
By definition,
\[ E(X)=\int_0^\infty x\lambda e^{-\lambda x}\,dx. \]
Use integration by parts with \(u=x\) and \(dv=\lambda e^{-\lambda x}dx\), so \(du=dx\) and \(v=-e^{-\lambda x}\). Then
\[ \begin{aligned} E(X) &=\left[-xe^{-\lambda x}\right]_0^\infty +\int_0^\infty e^{-\lambda x}\,dx\\ &=0+\frac1\lambda. \end{aligned} \]
lambda <- 0.5
mean_by_integration <- integrate(
function(x) x * dexp(x, rate = lambda),
lower = 0,
upper = Inf
)$value
c(integration = mean_by_integration,
formula = 1 / lambda)integration formula
2 2 A second integration by parts gives
\[ E(X^2)=\int_0^\infty x^2\lambda e^{-\lambda x}\,dx =2\int_0^\infty xe^{-\lambda x}\,dx =\frac{2}{\lambda^2}. \]
Therefore
\[ \operatorname{Var}(X) =E(X^2)-\{E(X)\}^2 =\frac{2}{\lambda^2}-\frac{1}{\lambda^2} =\frac{1}{\lambda^2}, \]
and
\[ \operatorname{SD}(X)=\frac1\lambda. \]
second_moment_by_integration <- integrate(
function(x) x^2 * dexp(x, rate = lambda),
lower = 0,
upper = Inf
)$value
variance_by_moments <-
second_moment_by_integration - mean_by_integration^2
c(
second_moment_integration = second_moment_by_integration,
second_moment_formula = 2 / lambda^2,
variance = variance_by_moments,
variance_formula = 1 / lambda^2,
standard_deviation = sqrt(variance_by_moments)
)second_moment_integration second_moment_formula variance
8 8 4
variance_formula standard_deviation
4 2 The CDF is
\[ F(x)=1-e^{-\lambda x}, \qquad x\geq0. \]
Solving \(F(Q_p)=p\) gives
\[ Q_p=-\frac{\log(1-p)}{\lambda}. \]
Thus
\[ Q_1=-\frac{\log(0.75)}{\lambda}, \qquad Q_2=-\frac{\log(0.50)}{\lambda}, \qquad Q_3=-\frac{\log(0.25)}{\lambda}. \]
For \(\lambda=0.5\), \(Q_1\approx0.5754\).
p <- c(Q1 = 0.25, Q2 = 0.50, Q3 = 0.75)
quantiles_formula <- -log(1 - p) / lambda
quantiles_R <- qexp(p, rate = lambda)
cbind(
probability = p,
formula = quantiles_formula,
R = quantiles_R
) probability formula R
Q1 0.25 0.5753641 0.5753641
Q2 0.50 1.3862944 1.3862944
Q3 0.75 2.7725887 2.7725887IQR_exponential <- unname(quantiles_R["Q3"] - quantiles_R["Q1"])
IQR_exponential[1] 2.197225