Random Variables — Turning Uncertainty Into Numbers

It is two in the morning and you are on call. Something is wrong with the API. You have a log file with 40,000 lines in it, one per request, each marked success or failure. Your manager sends one message: how bad is it?

Nobody wants 40,000 lines. What they want is a number — and then, quietly, a second number saying how much the first one moves around. "About one failure a minute" is useful. "About one failure a minute, but a bad minute has five" is far better, because it tells whoever is deciding whether to page the database team what the worst case looks like.

The move that gets you from a pile of outcomes to those two numbers is so small it looks like nothing: attach a number to every possible outcome, then reason about the numbers instead of the outcomes. That is a random variable. Everything else — averages, spread, correlation, the covariance matrix that principal component analysis pulls apart — is arithmetic on top of it.

A random variable is a function, not a variable

The name is misleading and it costs people weeks. A random variable is not a variable that happens to be random. It is a function: it takes an outcome of a random process and returns a number.

Write the process as a set of possible outcomes, the sample space $\Omega$. Flip three coins and $\Omega$ contains eight things: HHH, HHT, HTH, and so on. Those are not numbers, they are labels. A random variable $X$ is a rule that turns each label into a number:

Define $X$ = "number of heads" and you get $X(\text{HHH}) = 3$, $X(\text{HTH}) = 2$, $X(\text{TTT}) = 0$. Nothing random has happened yet. $X$ is a fixed, deterministic rule. The randomness lives in which outcome occurs; $X$ just translates.

A random variable is not a variable that is random. It is a fixed rule that assigns a number to every outcome, so that you can do arithmetic on uncertainty.

One process supports many random variables at once, and picking one is a modelling decision. From the same request log you could define $X$ = failures per minute, $Y$ = the slowest latency in the minute, or $Z$ = 1 if any failure occurred. Three functions, three different analyses.

By convention a capital $X$ is the rule and a lower-case $x$ is a value it took, so $P(X = 3)$ reads "the probability the rule returns 3".

Discrete and continuous

Random variables split into two kinds, and almost every mistake in this area comes from applying one kind's rules to the other.

• Discrete: the possible values can be listed, even if the list is infinite. Failures per minute (0, 1, 2, …), a class label, the number of tokens in a sentence, a die roll.
• Continuous: the possible values fill an interval. Latency in milliseconds, a temperature, a model's predicted probability, the weight on a neural network edge.

The dividing question is not "is it a whole number?" but "between any two possible values, is there always another possible value?" For latency, yes. For failure counts, no.

Discrete variables: the probability mass function

For a discrete variable, you can just write down the probability of each value. That list is the probability mass function, or PMF, written $p(x) = P(X = x)$. Here is one for failures per minute, measured from a week of logs:

A PMF has to obey exactly two rules, and they are both obvious once stated:

1. Nothing is negative: $p(x) \ge 0$ for every $x$. There is no such thing as a $-0.2$ chance.
2. It all adds to one: $\sum_x p(x) = 1$. Something must happen.

Check the table: $0.50 + 0.25 + 0.15 + 0.07 + 0.03 = 1.00$. If your probabilities do not sum to one you have a bug, not a distribution — assert it in code.

With the PMF you answer questions by adding: the chance of a minute with at least two failures is $0.15 + 0.07 + 0.03 = 0.25$, one minute in four.

Continuous variables: density is not probability

Now try the same thing for latency. What is the probability a request takes exactly 137.0000… milliseconds? Zero — not "very small", exactly zero. Any interval holds infinitely many real numbers, so no single one can carry probability without the total blowing past one.

So continuous variables get a different object: the probability density function, or PDF, written $f(x)$. Probability is not the height of $f$; it is the area under $f$:

with $f(x) \ge 0$ everywhere and $\int_{-\infty}^{\infty} f(x)\,dx = 1$ — the same two axioms, with the sum replaced by an integral.

For a continuous variable, $f(x)$ is not a probability. Probability is area under $f$, and the area over a single point is zero. This is the single most common misunderstanding in the whole subject.

Here is the demonstration that usually fixes it for good. If latency is uniform between 0 and 200 ms, the density is flat at $f(x) = 1/200 = 0.005$ per millisecond — comfortably under one. But if some quantity is uniform between 0 and 0.5, the density must integrate to one over a width of 0.5, so $f(x) = 2$. A density of 2. If densities were probabilities that would be nonsense. The units are "probability per unit of $x$", and squeezing the same total probability into a narrower interval makes the density taller.

For the 0–200 ms case, $P(100 \le X \le 150) = 50 \times 0.005 = 0.25$ — width times height, because the density is flat. An integral just does the same thing for a curve that is not flat.

The CDF answers every question you actually ask

The cumulative distribution function is one function that works identically for both kinds:

"What is the probability the value lands at $x$ or below?" Three properties fall straight out of that: it never decreases, it starts at 0 far to the left, and it ends at 1 far to the right.

Two things make it the workhorse. First, interval probabilities become a subtraction:

From the failures table, $P(1 \le X \le 3) = F(3) - F(0) = 0.97 - 0.50 = 0.47$. Cross-check by adding the PMF entries: $0.25 + 0.15 + 0.07 = 0.47$. They agree, as they must.

Second, running it backwards gives percentiles. If latency is uniform on $[0, 200]$ then $F(x) = x/200$, so the 95th percentile — the p95 on your monitoring dashboard — solves $x/200 = 0.95$, giving $x = 190$ ms. That inverse is the quantile function, spelled .ppf() in SciPy.

Expectation, and the surprising power of linearity

The expected value is the average value you would see if you could repeat the process forever, with each possible value weighted by how likely it is:

On the failures table:

About 0.88 failures per minute. Be honest about what that means: you will never observe 0.88 failures in a minute, because it is not a possible value. Expectation is a long-run average, a balance point — not a prediction of any single observation. Forget that and you promise stakeholders "the expected wait is 4 minutes", then get baffled when half the users wait longer.

Linearity holds even when variables are tangled together

Expectation obeys a rule that is far stronger than it looks:

The remarkable part is the middle term. It holds whether or not $X$ and $Y$ are independent. Nothing else here is that forgiving. You can chop a horribly tangled quantity into simple pieces, average each piece on its own, and add — without ever working out how the pieces interact.

Here is a case where that saves real work. Five items in a batch, and a model that outputs the five correct labels in a random order. Let $X$ be the number of positions where the shuffled prediction lands on the right label. How many matches do you expect?

Getting the distribution of $X$ directly means counting permutations with exactly $k$ fixed points, which is unpleasant. Linearity sidesteps it. Define an indicator $I_i$ equal to 1 if position $i$ matches and 0 otherwise, so $X = I_1 + \cdots + I_5$. Each position has a $1/5$ chance of receiving its own label, so $E[I_i] = 0.2$ and

One match expected, on average. The indicators are strongly dependent — if four positions match, the fifth is forced to match too — and linearity did not care. Grinding out all 120 permutations confirms it: 44 have no match, 45 have one, 20 have two, 10 have three and 1 has all five, which averages to exactly $120/120 = 1$.

Variance, standard deviation, and when they add

The average alone hides the thing you were actually worried about at 02:00: how bad does a bad minute get? Variance measures spread as the average squared distance from the mean:

Squaring does two jobs: it stops positive and negative deviations cancelling, and it punishes large deviations disproportionately. The second form is faster to compute. Both should agree — check, with $\mu = 0.88$.

By the definition, take each deviation, square it, weight it:

By the shortcut, first get $E[X^2] = 0(0.50) + 1(0.25) + 4(0.15) + 9(0.07) + 16(0.03) = 0.25 + 0.60 + 0.63 + 0.48 = 1.96$, then subtract:

Identical. Note that $E[X^2] = 1.96$ is not the same as $(E[X])^2 = 0.7744$ — the gap between them is the variance.

Variance is in squared units, and "1.19 failures squared" means nothing to a human. Take the square root to get the standard deviation, $\sigma = \sqrt{1.1856} \approx 1.09$ failures — back in the original units and directly comparable to the mean of 0.88. That is the number you quote: typical minute 0.88, typical wobble 1.09, so three or four failures is unremarkable and forty is an incident.

The two rules, and the one that has a condition

Two consequences worth keeping. Adding a constant does nothing to variance, since everything shifts together: $\mathrm{Var}(X + 100) = \mathrm{Var}(X)$. And $\mathrm{Var}(2X + 3) = 4 \times 1.1856 = 4.7424$ — doubling a quantity quadruples its variance.

Expectation always adds. Variance only adds when the variables are independent — otherwise you owe a covariance term, and forgetting it is how error bars end up far too narrow.

Two variables at once: joint, marginal, covariance, correlation

Real data has many columns. The joint distribution $P(X = x, Y = y)$ gives the probability of both values together. Sum over one variable and you recover a marginal — the distribution of the other on its own: $P(X = x) = \sum_y P(X = x, Y = y)$.

$X$ and $Y$ are independent when knowing one tells you nothing about the other. The test is crisp: independence holds if and only if $P(X, Y) = P(X)P(Y)$ for every pair.

Covariance

Covariance asks whether two variables move together:

When both are above their means, or both below, the product is positive; when one is up and the other down, it is negative. Averaging those products gives the direction of the relationship. Take five training runs — hours of training $x$ against validation accuracy $y$:

Deviations: $x - \bar{x} = [-3, -1, -1, 1, 4]$ and $y - \bar{y} = [-14, -6, -2, 4, 18]$. Products: $42, 6, 2, 4, 72$, summing to 126. Dividing by $n = 5$:

Positive — more training goes with higher accuracy. But 25.2 of what? The units are hours times accuracy-points. Measure training in minutes instead and every $x$ deviation multiplies by 60, so the covariance becomes 1512. The relationship did not change, only the ruler. Covariance gives you the sign reliably and the strength not at all.

Correlation fixes the units

Divide by both standard deviations and the units cancel:

Here $\sigma_x = \sqrt{28/5} = \sqrt{5.6} \approx 2.366$ and $\sigma_y = \sqrt{576/5} = \sqrt{115.2} \approx 10.733$, so

Correlation is always between $-1$ and $+1$, and it is unitless: switch hours to minutes and it stays 0.992. That is why every report quotes correlation, never covariance.

The trap: $\rho = 0$ does not mean independent

Correlation measures linear association only. Let $X$ be uniform on $\{-2, -1, 0, 1, 2\}$ and let $Y = X^2$ exactly — so $Y$ is $\{4, 1, 0, 1, 4\}$, completely determined by $X$, as dependent as two variables can possibly be. Compute the covariance. $\mu_X = 0$ and $\mu_Y = (4+1+0+1+4)/5 = 2$, so

Correlation zero, dependence total — the two arms of the parabola cancel exactly.

Zero correlation means no linear relationship. It does not mean no relationship. Dropping a feature because its correlation with the target is near zero can throw away the most predictive column you have.

The covariance matrix

With $d$ features, collect every pairwise covariance into a $d \times d$ matrix $\Sigma$, where $\Sigma_{ij} = \mathrm{Cov}(X_i, X_j)$. The diagonal holds the variances — how much each feature varies on its own, since $\mathrm{Cov}(X_i, X_i) = \mathrm{Var}(X_i)$. The off-diagonal entries hold the redundancy — how much each pair moves together. It is symmetric, because covariance does not care about order.

This is the object principal component analysis takes apart: its eigenvectors are the directions along which the data varies most, its eigenvalues say how much variance sits along each. When you hear "PCA finds the directions of maximum variance", $\Sigma$ is what gets decomposed.

Doing all of this in NumPy

1import numpy as np23# --- discrete variable: failures per minute ---4values = np.array([0, 1, 2, 3, 4])5probs  = np.array([0.50, 0.25, 0.15, 0.07, 0.03])6assert np.isclose(probs.sum(), 1.0)      # axiom 2, as a runtime check78mu   = (values * probs).sum()            # 0.889ex2  = (values**2 * probs).sum()         # 1.9610var  = ex2 - mu**2                       # 1.185611print(mu, var, np.sqrt(var))             # 0.88 1.1856 1.0888...1213cdf = np.cumsum(probs)                   # [0.5 0.75 0.9 0.97 1.0]14print(cdf[3] - cdf[0])                   # P(1 to 3 inclusive) = 0.47

1# --- covariance and correlation on paired data ---2hours = np.array([2, 4, 4, 6, 9])3acc   = np.array([60, 68, 72, 78, 92])45# np.cov defaults to ddof=1 (the sample estimator, dividing by n-1).6# Pass ddof=0 when you want the population figure worked above.7print(np.cov(hours, acc, ddof=0))8# [[  5.6  25.2]9#  [ 25.2 115.2]]   diagonal = variances, off-diagonal = covariance1011print(np.corrcoef(hours, acc)[0, 1])           # 0.99221213print(np.cov(hours * 60, acc, ddof=0)[0, 1])   # 1512.0  -- covariance moved14print(np.corrcoef(hours * 60, acc)[0, 1])      # 0.9922  -- correlation did not

1# --- the zero-correlation trap, and a full covariance matrix ---2x = np.array([-2, -1, 0, 1, 2])3y = x ** 24print(np.corrcoef(x, y)[0, 1])           # 0.0, yet y is a function of x56rng = np.random.default_rng(0)7X = rng.normal(size=(1000, 3)) * np.array([1.0, 0.5, 3.0])8Sigma = np.cov(X, rowvar=False)          # rowvar=False: rows are samples9print(np.diag(Sigma).round(2))           # ~[1.0 0.25 9.0] variances10print(np.sqrt(np.diag(Sigma)).round(2))  # ~[1.0 0.5  3.0] std deviations1112# --- Monte Carlo: get the same numbers by simulating instead ---13draws = rng.choice(values, size=200_000, p=probs)14print(draws.mean(), draws.var())         # ~0.88, ~1.186

That last trick matters. When a distribution is too awkward to integrate, sample from it a few hundred thousand times and average. That is Monte Carlo estimation, and it is how dropout uncertainty, reinforcement-learning returns and Bayesian posteriors get evaluated in practice.

What this changes when you build a model

Every one of these ideas shows up as a line of code in a real pipeline, and knowing which one you are using stops a specific mistake.

Feature scaling is expectation and variance. Standardisation computes $z = (x - \mu)/\sigma$ per column, and the rule $\mathrm{Var}(aX + b) = a^2\mathrm{Var}(X)$ is exactly why that leaves every feature with variance 1. It also explains why $\mu$ and $\sigma$ must come from the training set and be reused unchanged on the test set: they are parameters estimated from data, and recomputing them on test data leaks information.

The covariance matrix is a diagnostic, not just PCA input. Print $\Sigma$ before you train anything. Large off-diagonal entries mean redundant features, which makes linear-model coefficients unstable and uninterpretable. A near-zero diagonal entry is a constant column that contributes nothing and breaks anything dividing by $\sigma$.

Your model's accuracy is a random variable. Retrain with a different seed, shuffle or split and the number changes. Reporting "87.3%" with no sense of its spread quotes a single draw as if it were a constant — which is how teams ship a model that is not actually better than the one it replaced.

Bias and variance are literally these quantities. Treat the prediction at a point as a random variable across possible training sets. Bias is how far its expectation sits from the truth; variance is how much it moves as the training set changes. A deep tree has low bias and high variance; a linear model usually has the opposite. Every regularisation technique trades one against the other.

So the next time someone asks "how bad is it?", you have a real answer: define the random variable, compute its expectation and its standard deviation, and quote both. The mean tells people what to plan for; the spread tells them what to survive.