---
title: "Expectation, Variance, and Quantiles"
subtitle: "SD04 — Random Variables"
author: "Samer El Kababji"
output:
  html_document:
    toc: true
format:
  html:
    code-links:
      - text: "Download R Markdown"
        icon: file-earmark-arrow-down
        href: SD04-02-expectation-variance-and-quantiles.Rmd?download=1
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = TRUE)
```

**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.

## Expected value

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).
$$

```{r discrete-expectation-function}
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)
)
```

For a continuous random variable $X$ with PDF $f$,

$$
E(X)=\mu_X=\int_{-\infty}^{\infty}x f(x)\,dx.
$$

```{r continuous-expectation-function}
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)
```

The 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 and standard deviation

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)}.
$$

```{r discrete-moment-function}
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))
```

## Percentiles and quantiles

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.

```{r discrete-quantile-function}
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)
)
```

## Worked example: machine-breakdown cost

Let $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.
$$

```{r machine-breakdown-event}
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
```

The 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.
$$

```{r machine-breakdown-moments}
breakdown_moments <- discrete_moments(
  breakdown_cost,
  breakdown_probability
)

breakdown_moments
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.
$$

```{r machine-breakdown-quantiles}
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"]
)
```

These are **distribution quantiles**, obtained from the PMF. They should not
be confused with an interpolated sample-quantile convention applied to a data
vector.

## Worked example: exponential random variable

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}
$$

### Mean

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}
$$

```{r exponential-mean}
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)
```

### Second moment, variance, and standard deviation

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.
$$

```{r exponential-moments}
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)
)
```

### Exponential quantiles

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$.

```{r exponential-quantiles}
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
)

IQR_exponential <- unname(quantiles_R["Q3"] - quantiles_R["Q1"])
IQR_exponential
```

## Practice questions

1. A discrete random variable takes values $-1,0,2$ with probabilities
   $0.2,0.5,0.3$. Find its mean, variance, and standard deviation.
2. For the machine-breakdown distribution, find $E(2X+10)$ directly from
   its PMF.
3. Explain why the machine-breakdown median is 200 even though exactly 20% of
   the probability mass occurs at 200.
4. For $X\sim\operatorname{Exp}(0.25)$, find its mean, variance, and median.
5. Show that the exponential IQR is $\log(3)/\lambda$.

::: {.callout-note collapse="true"}
### Check your answers

1. $E(X)=0.4$, $E(X^2)=1.4$, so
   $\operatorname{Var}(X)=1.4-0.4^2=1.24$ and
   $\operatorname{SD}(X)\approx1.114$.
2. $E(2X+10)=2(230)+10=470$; direct PMF weighting gives the same result.
3. The median is the first support value whose CDF reaches 0.5:
   $F(50)=0.3<0.5$ and $F(200)=0.5$.
4. The mean is 4, the variance is 16, and the median is
   $-\log(0.5)/0.25\approx2.773$.
5. $Q_3-Q_1=\{-\log(0.25)+\log(0.75)\}/\lambda
   =\log(3)/\lambda$.
:::