Expectation, Variance, and the CLT

The expected value (mean) of $X$ is its probability-weighted average:

$\mathbb{E}[X] = \begin{cases} \displaystyle\sum_x x\, p(x) & \text{discrete} \\ \displaystyle\int_{-\infty}^{\infty} x\, f(x)\, dx & \text{continuous} \end{cases}$

Think of it as the center of mass of the PMF/PDF. It need not be a likely value — for a fair die, $\mathbb{E}[X] = 3.5$, which never appears on a roll.

Expected value of a function: $\mathbb{E}[g(X)] = \sum_x g(x)\, p(x)$. This is the law of the unconscious statistician (LOTUS) — you integrate $g(x)$ against the original distribution of $X$, no change of variables needed.

The single most useful fact in probability:

$\mathbb{E}[aX + bY + c] = a\,\mathbb{E}[X] + b\,\mathbb{E}[Y] + c$

This holds regardless of whether $X$ and $Y$ are independent. No joint distribution needed.

Example: what is the expected number of fixed points (cards in their natural position) in a random shuffle of $n$ cards? Define $X_i = \mathbf{1}[\text{card } i \text{ is fixed}]$. Then $\mathbb{E}[X_i] = 1/n$ and $\mathbb{E}\!\left[\sum_{i=1}^n X_i\right] = \sum_{i=1}^n \mathbb{E}[X_i] = n \cdot \frac{1}{n} = 1.$ Always exactly 1, for any $n$ — a clean result that would be a combinatorial nightmare without linearity.

Variance measures spread around the mean:

$\text{Var}(X) = \mathbb{E}\!\left[(X - \mathbb{E}[X])^2\right] = \mathbb{E}[X^2] - (\mathbb{E}[X])^2.$

The second form — raw second moment minus squared mean — is almost always faster to compute. The standard deviation is $\sigma = \sqrt{\text{Var}(X)}$, which lives in the same units as $X$.

Key scaling rules:

• $\text{Var}(aX + b) = a^2\,\text{Var}(X)$ — shift has no effect, scale is squared.
• If $X, Y$ independent: $\text{Var}(X + Y) = \text{Var}(X) + \text{Var}(Y)$.
• If dependent: $\text{Var}(X + Y) = \text{Var}(X) + \text{Var}(Y) + 2\,\text{Cov}(X, Y)$.

Variance is not linear — a common source of errors in error-propagation and mini-batch gradient analysis.

Let $X$ and $Y$ be independent fair dice ($1$–$6$). Compute $\mathbb{E}[X + Y]$ and $\text{Var}(X + Y)$.

Step 1. For one fair die: $\mathbb{E}[X] = \frac{1}{6}(1+2+3+4+5+6) = 3.5$.

Step 2. $\mathbb{E}[X^2] = \frac{1}{6}(1+4+9+16+25+36) = \frac{91}{6} \approx 15.17$.

Step 3. $\text{Var}(X) = 15.17 - 3.5^2 = 15.17 - 12.25 = 2.92$.

Step 4. By linearity: $\mathbb{E}[X + Y] = 3.5 + 3.5 = 7$.

Step 5. By independence: $\text{Var}(X + Y) = 2.92 + 2.92 = 5.83$, so $\sigma_{X+Y} = \sqrt{5.83} \approx 2.42$.

import numpy as np
rolls = np.random.randint(1, 7, size=(10_000_000, 2))
sums = rolls.sum(axis=1)
print(f"E[X+Y] = {sums.mean():.4f}")  # ~7.0000
print(f"Var    = {sums.var():.4f}")   # ~5.8333

Covariance quantifies how two RVs move together:

$\text{Cov}(X, Y) = \mathbb{E}[(X - \mathbb{E}[X])(Y - \mathbb{E}[Y])] = \mathbb{E}[XY] - \mathbb{E}[X]\,\mathbb{E}[Y].$

$\text{Cov}(X, X) = \text{Var}(X)$. Covariance is scale-dependent, making it hard to interpret directly.

Pearson correlation normalizes it: $\rho(X, Y) = \frac{\text{Cov}(X, Y)}{\sqrt{\text{Var}(X)\,\text{Var}(Y)}} \in [-1, 1].$

$\rho = 1$: perfect positive linear relationship; $\rho = 0$: uncorrelated (but not necessarily independent); $\rho = -1$: perfect negative.

For a random vector $\mathbf{X} \in \mathbb{R}^d$, the covariance matrix $\Sigma_{ij} = \text{Cov}(X_i, X_j)$ is symmetric and positive semi-definite — it appears in Gaussian distributions, PCA, and Kalman filters.

Let $X_1, X_2, \ldots$ be i.i.d. with mean $\mu$ and finite variance $\sigma^2$. The sample mean is $\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i$.

Weak LLN: For any $\varepsilon > 0$, $P(|\bar{X}_n - \mu| > \varepsilon) \to 0$ as $n \to \infty$.

Proof sketch via Chebyshev: $P(|\bar{X}_n - \mu| > \varepsilon) \le \frac{\text{Var}(\bar{X}_n)}{\varepsilon^2} = \frac{\sigma^2}{n\varepsilon^2} \to 0$.

In ML: Monte Carlo estimators, stochastic gradient descent, empirical risk minimization — all implicitly rely on the LLN. When you average enough noisy gradient estimates, you get close to the true gradient. The LLN tells you convergence happens; the CLT tells you how fast and what the fluctuations look like.

The most important theorem in statistics:

CLT. If $X_1, \ldots, X_n$ are i.i.d. with mean $\mu$ and variance $\sigma^2 < \infty$, then $\frac{\bar{X}_n - \mu}{\sigma / \sqrt{n}} \xrightarrow{d} \mathcal{N}(0, 1) \quad \text{as } n \to \infty.$

Three takeaways:

1. The shape of the original distribution is irrelevant — the sum converges to Gaussian.
2. The fluctuations shrink as $1/\sqrt{n}$, not $1/n$. Halving the error requires 4× the data.
3. The approximation is often usable by $n \approx 30$ for unimodal distributions, but may need $n \gg 100$ for heavy-tailed ones.

This is why Gaussian noise appears in so many models: real noise is typically a sum of many small independent perturbations.

The standard error of the mean (SEM) is the standard deviation of $\bar{X}_n$: $\text{SE} = \frac{\sigma}{\sqrt{n}}.$

In practice $\sigma$ is unknown, so you substitute the sample standard deviation $s = \sqrt{\frac{1}{n-1}\sum_i (X_i - \bar{X})^2}$.

The SEM quantifies estimation uncertainty. A rough 95% interval for $\mu$ is $\bar{X}_n \pm 1.96 \cdot \text{SE}$ (via the CLT). Comparing two models? Compute the SEM of the accuracy difference — if the interval excludes 0, the difference is statistically detectable at the 5% level.

import numpy as np
from scipy import stats

np.random.seed(42)
data = np.random.exponential(scale=5, size=200)  # non-Gaussian!
xbar = data.mean()
se   = stats.sem(data)
ci   = stats.t.interval(0.95, df=len(data)-1, loc=xbar, scale=se)
print(f"mean={xbar:.2f}, SE={se:.2f}, 95% CI={ci}")

Notice we used a $t$-distribution (not normal) because $n=200$ with unknown $\sigma$ — conservative and correct.

How repeated averaging turns any distribution into a Gaussian:

Beyond the CLT, Gaussians dominate ML for three deeper reasons:

1. Maximum entropy. Among all distributions with fixed mean $\mu$ and variance $\sigma^2$, $\mathcal{N}(\mu, \sigma^2)$ has the highest differential entropy. It's the "most uncertain" distribution given those two constraints — a natural prior when you know only the first two moments.

2. Analytic tractability. Gaussians are closed under linear transformations, conditioning, and marginalization. If $(X, Y)$ is jointly Gaussian, $X \mid Y=y$ is Gaussian with a closed-form mean and variance. This is the engine of Gaussian processes and Kalman filters.

3. Exponential family. $\mathcal{N}(\mu, \sigma^2)$ is in the exponential family, making maximum-likelihood estimation and Bayesian conjugacy computationally clean.

Every time you assume "Gaussian noise," you're implicitly invoking the CLT, maximum entropy, and analytic convenience simultaneously.