1. Foundations I: Random Variables, Expectations, Quantiles, and Convexity

Learning Objectives. By the end of this chapter, you will be able to translate uncertain physical events into Random Variables ($X: \Omega \to \mathbb{R}$) and Empirical Cumulative Distribution Functions ($\hat{F}_n(x)$); compute and contrast Arithmetic Expectation ($\mathbb{E}[X]$) against Order Statistics and Quantiles ($q_{0.05}, q_{0.50}, q_{0.95}$); interpret Variance ($\sigma^2$), Covariance ($\mathrm{Cov}(X,Y)$), and Pearson Correlation ($\rho_{X,Y}$) geometrically as inner products of centered vectors; derive Convexity ($f’’(x) \ge 0$) and Jensen’s Inequality ($\mathbb{E}[f(X)] \ge f(\mathbb{E}[X])$) from first principles; and distinguish vertical distribution distance (Continuous Ranked Probability Score, $\mathrm{CRPS}$) from horizontal transport distance (Wasserstein-1, $W_1$).

Step 1 — Intuition and the Real-World Paradox

The River With an Average Depth of One Meter

A civil engineer is asked whether an infantry platoon—whose soldiers have an average height of $1.75\text{ m}$ and carry radio packs strapped at $1.40\text{ m}$—can safely wade across a river whose hydrographic bulletin reports:

“Mean depth across the 100-meter crossing: $\bar{x} = 1.00\text{ m}$.”

Relying on the reported mean, the engineer approves the crossing. Halfway across, the riverbed drops into a 20-meter navigational channel that is $2.40\text{ m}$ deep. Every radio is flooded and the platoon is swept downstream.

Why did a true mathematical statement ($\bar{x} = 1.00\text{ m}$) produce an operational catastrophe? Because two structural operations happened silently between the river and the decision:

1. Collapsing a Distribution into a Single Number: The river’s depth is not a scalar; it is a distribution across space ($80\text{ m}$ of $0.65\text{ m}$ shallows and $20\text{ m}$ of $2.40\text{ m}$ channel: $0.80 \times 0.65 + 0.20 \times 2.40 = 1.00\text{ m}$).
2. Evaluating a Nonlinear Payoff at the Mean: The operational damage function $f(x)$—radio destruction and drowning—is not linear in depth $x$. Going from $0.30\text{ m}$ to $0.60\text{ m}$ causes zero damage ($f(0.65) = 0$), whereas crossing $1.40\text{ m}$ causes $100%$ equipment loss ($f(2.40) = 100$).

Notice the difference between the two numbers:

• Damage at the average depth: $f(\mathbb{E}[X]) = f(1.00\text{ m}) = 0%$ loss.
• Average damage across the actual river: $\mathbb{E}[f(X)] = 0.80 \times f(0.65) + 0.20 \times f(2.40) = 0.80(0) + 0.20(100) = 20%$ loss (and $100%$ failure for anyone walking the full crossing!).

Before we can understand the Three Ontological Layers of Chapter 2 ($L_1 \to L_2 \to L_3$)—including Meng’s correlation identity, Jensen’s Gap, $\mathrm{CRPS}$, and Wasserstein-1 performative drift—we must build the four elementary pillars of quantitative stochology from scratch:

1. Random Variables, CDFs, and Quantiles (how a full distribution is represented in memory);
2. Moments, Variance, Covariance, and Correlation (how dispersion and co-movement are measured);
3. Convexity and Secant Geometry (why curves bend averages);
4. Distances Between Distributions ($\mathrm{CRPS}$ vs. Wasserstein-1 $W_1$ and Lipschitz sensitivity).

Step 2 — Foundations and Mathematical Formalization

1.1 From Outcomes to Random Variables, CDFs, and Quantiles

In the Preface, we worked with discrete events $A \subseteq \Omega$ (“the storm hits,” “the test is positive”). In engineering and science, however, observations are numerical measurements: wind speed in knots, queue wait time in seconds, financial loss in millions of krone.

Definition 1.1 (Random Variable and Cumulative Distribution Function)

A random variable $X: \Omega \to \mathbb{R}$ is a deterministic rule that assigns a real number $X(\omega)$ to every elementary outcome $\omega \in \Omega$. Its complete probabilistic law is captured by its Cumulative Distribution Function (CDF) $F_X: \mathbb{R} \to [0, 1]$:

$$F_X(x) = P({\omega \in \Omega : X(\omega) \le x}) = P(X \le x)$$

Every CDF satisfies three properties:

1. Monotonicity: If $a \le b$, then $F_X(a) \le F_X(b)$ (adding more possible values can never shrink probability).
2. Boundary limits: $\lim_{x \to -\infty} F_X(x) = 0$ and $\lim_{x \to +\infty} F_X(x) = 1$.
3. Interval probability: For any interval $(a, b]$, $P(a < X \le b) = F_X(b) - F_X(a)$.

Empirical CDF and Order Statistics (QuantileSummary)

When we observe a finite sample of $n$ measurements $x_1, x_2, \dots, x_n$ in Rust, we first sort them in ascending order to obtain the order statistics:

$$x_{(1)} \le x_{(2)} \le \dots \le x_{(n)}$$

The Empirical CDF $\hat{F}_n(x)$ is the staircase function that jumps by $1/n$ at each sorted data point:

$$\hat{F}_n(x) = \frac{1}{n} \sum_{i=1}^n \mathbf{1}{x_i \le x}$$

Inverting $\hat{F}_n$ gives the $\alpha$-quantile $q_\alpha = F_X^{-1}(\alpha) = \inf{x \in \mathbb{R} : F_X(x) \ge \alpha}$ for any probability level $\alpha \in (0, 1)$:

• $q_{0.05}$ (5th percentile): Only $5%$ of observations fall at or below this value (lower tail).
• $q_{0.50}$ (Median): Exactly half the probability mass lies below and half above.
• $q_{0.95}$ (95th percentile): $95%$ of observations fall below, and $5%$ exceed this threshold (upper tail).

Worked Example 1.1: Why the Median Ignores Catastrophic Outliers while the Mean Explodes

Consider five daily server latency readings (in ms):

$$x = [10, ; 12, ; 14, ; 16, ; 18]$$

• The sample mean is $\bar{x} = \frac{10 + 12 + 14 + 16 + 18}{5} = 14.0\text{ ms}$.
• The sample median ($q_{0.50}$, middle element $x_{(3)}$) is $14.0\text{ ms}$.

Now suppose on Friday the server suffers a memory stall, so the fifth reading is $1\text{,}018\text{ ms}$ instead of $18\text{ ms}$:

$$x’ = [10, ; 12, ; 14, ; 16, ; 1\text{,}018]$$

• The median $x’_{(3)}$ is still $14.0\text{ ms}$—it depends only on rank order, so a single extreme observation cannot move it past its neighbors.
• The mean jumps to $\bar{x}’ = \frac{1\text{,}070}{5} = \mathbf{214.0\text{ ms}}$, which describes neither the four normal days ($10$–$16\text{ ms}$) nor the stalled day ($1\text{,}018\text{ ms}$)!

1.2 Expectation, Variance, Covariance, and Pearson Correlation

Expectation as a Weighted Sum (Linear Operator)

For a discrete random variable taking values $x_1, \dots, x_k$ with probabilities $p_i = P(X = x_i)$ (or a continuous random variable with density $p(x) = F’_X(x)$), the expectation (or expected value) is the probability-weighted center of mass:

$$\mathbb{E}[X] = \sum_{i=1}^k x_i , p_i \qquad \text{or} \qquad \mathbb{E}[X] = \int_{-\infty}^{\infty} x , p(x),dx$$

The single most important algebraic property of $\mathbb{E}[\cdot]$ is linearity: for any constants $a, b \in \mathbb{R}$ and any two random variables $X, Y$ (even if they are strongly dependent!),

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

Variance and Standard Deviation

How spread out is $X$ around its center $\mu_X = \mathbb{E}[X]$? If we averaged the raw deviations $\mathbb{E}[X - \mu_X]$, positive and negative deviations would cancel to zero ($\mathbb{E}[X] - \mu_X = 0$). Squaring the deviations before averaging gives the variance:

$$\mathrm{Var}(X) = \sigma_X^2 = \mathbb{E}!\left[(X - \mu_X)^2\right] = \mathbb{E}[X^2] - (\mathbb{E}[X])^2$$

Because $\mathrm{Var}(X)$ has squared units (e.g., $\text{ms}^2$ or $\text{dollars}^2$), taking its square root restores the original physical units, giving the standard deviation $\sigma_X = \sqrt{\mathrm{Var}(X)}$.

Covariance and Pearson Correlation as Vector Cosine ($\rho_{X,Y}$)

Suppose we observe two variables simultaneously on each unit $i \in {1, \dots, n}$: for example, $X_i$ = wind speed and $Y_i$ = storm surge, or (as we will see in Chapter 2 for Meng’s Identity) $R_i \in {0, 1}$ = whether unit $i$ was recorded in our database and $Y_i$ = its true outcome.

Define the centered vectors in $\mathbb{R}^n$ by subtracting their sample means:

$$\tilde{x}_i = x_i - \bar{x}, \qquad \tilde{y}_i = y_i - \bar{y}$$

1. The sample covariance is the average dot product of these centered vectors: $$\mathrm{Cov}(X, Y) = \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y}) = \frac{1}{n}\langle \tilde{\mathbf{x}}, \tilde{\mathbf{y}} \rangle$$ When $x_i > \bar{x}$ and $y_i > \bar{y}$ happen together (or both fall below their means together), the product $(x_i - \bar{x})(y_i - \bar{y})$ is positive.
2. The Pearson correlation coefficient $\rho_{X,Y}$ normalizes covariance by the product of the two standard deviations, which is identically the cosine of the angle $\theta$ between the centered vectors $\tilde{\mathbf{x}}$ and $\tilde{\mathbf{y}}$ in $\mathbb{R}^n$:

$$\rho_{X,Y} = \frac{\mathrm{Cov}(X,Y)}{\sigma_X \sigma_Y} = \frac{\langle \tilde{\mathbf{x}}, \tilde{\mathbf{y}} \rangle}{|\tilde{\mathbf{x}}|_2 , |\tilde{\mathbf{y}}|_2} = \cos(\theta) \in [-1, 1]$$

By the Cauchy-Schwarz Inequality ($|\langle \tilde{\mathbf{x}}, \tilde{\mathbf{y}} \rangle| \le |\tilde{\mathbf{x}}|_2 |\tilde{\mathbf{y}}|_2$), correlation is strictly bounded between $-1$ and $+1$. In Chapter 2, this simple inner-product identity will reveal why non-random missing data ($R_i$ correlated with $Y_i$, $\rho_{R,Y} \ne 0$) destroys Big Data averages!

1.3 Nonlinear Functions, Convexity, and the Derivation of Jensen’s Inequality

Why does $\mathbb{E}[f(X)]$ almost never equal $f(\mathbb{E}[X])$ when $f$ is curved?

Definition 1.2 (Convex Function)

A function $f: \mathbb{R} \to \mathbb{R}$ is convex if, for every pair of points $x_1, x_2 \in \mathbb{R}$ and every weight $\lambda \in [0, 1]$, the straight chord (secant line) connecting $(x_1, f(x_1))$ and $(x_2, f(x_2))$ lies on or above the graph of $f$:

$$\lambda f(x_1) + (1 - \lambda) f(x_2) \ge f!\left(\lambda x_1 + (1 - \lambda) x_2\right)$$

If $f$ is twice differentiable, $f$ is convex if and only if its curvature (second derivative) is non-negative everywhere:

$$f’’(x) \ge 0 \quad \forall x \in \mathbb{R}$$

(Conversely, if $f’’(x) \le 0$, $f$ is concave, like $\ln(x)$ or $\sqrt{x}$, and the inequality reverses: $\mathbb{E}[f(X)] \le f(\mathbb{E}[X])$.)

Step-by-Step Proof of Jensen’s Inequality via Supporting Tangent Lines

Let $\mu = \mathbb{E}[X]$. Because a convex function $f$ always curves upward away from its tangent line at $x = \mu$, write the second-order Taylor expansion of $f(X)$ around $\mu$:

$$f(X) = f(\mu) + f’(\mu)(X - \mu) + \frac{1}{2} f’’(\xi)(X - \mu)^2$$

for some intermediate point $\xi$ between $\mu$ and $X$. Now take the expectation $\mathbb{E}[\cdot]$ of both sides using linearity:

$$\mathbb{E}[f(X)] = f(\mu) + f’(\mu)\underbrace{\mathbb{E}[X - \mu]}_{= , 0} + \frac{1}{2}\mathbb{E}!\left[f’’(\xi)(X - \mu)^2\right]$$

Look at the three terms on the right-hand side:

1. $f(\mu) = f(\mathbb{E}[X])$ is the naïve evaluation at the mean.
2. The first-order linear slope term $f’(\mu)\mathbb{E}[X - \mu]$ vanishes identically because $\mathbb{E}[X - \mu] = \mu - \mu = 0$.
3. Because $f’’(\xi) \ge 0$ and $(X - \mu)^2 \ge 0$, the remaining curvature term is non-negative! Calling the difference Jensen’s Gap ($\Delta_J$), we obtain:

$$\Delta_J = \mathbb{E}[f(X)] - f(\mathbb{E}[X]) \approx \frac{1}{2} f’’(\mu),\mathrm{Var}(X) \ge 0$$

The Fundamental Lesson of Jensen’s Gap: Jensen’s Gap is proportional to curvature times variance ($\frac{1}{2} f’‘(\mu)\sigma_X^2$). Whenever your operational cost function has positive curvature ($f’’ > 0$, such as hospital queue delays, option payoffs, or hurricane damage $f(v) \propto v^3$) OR your system has large variance ($\sigma_X^2 \gg 0$), evaluating your formula at the average input $f(\mathbb{E}[X])$ systematically understates the true average cost $\mathbb{E}[f(X)]$.

1.4 Comparing Full Distributions: CRPS vs. Wasserstein-1 ($W_1$) and Lipschitz Continuity

Once we refuse to collapse a random variable $X$ into a lone mean, how do we measure the error of a probabilistic forecast $F(x)$ or the distance between two distributions $P$ and $Q$?

1. Vertical Distance: Continuous Ranked Probability Score ($\mathrm{CRPS}$)

Suppose a forecaster issues a full predictive CDF $F(x) = P(X \le x)$, and nature subsequently reveals the actual scalar outcome $y \in \mathbb{R}$. Nature’s outcome can itself be written as a step-function CDF $\mathbf{1}{x \ge y}$ that jumps from $0$ to $1$ at $x = y$.

The Continuous Ranked Probability Score (crps_empirical in crates/stochology-core/src/layers.rs) integrates the squared vertical gap between the forecast CDF $F(x)$ and nature’s step function across the entire real line:

$$\mathrm{CRPS}(F, y) = \int_{-\infty}^{\infty} \left(F(x) - \mathbf{1}{x \ge y}\right)^2 dx = \mathbb{E}_F[|X - y|] - \frac{1}{2}\mathbb{E}_F[|X - X’|]$$

where $X, X’$ are independent draws from $F$. Notice the two parts of the right-hand representation:

• $\mathbb{E}_F[|X - y|]$ penalizes the forecast if its samples are far from the true outcome $y$ (accuracy).
• $-\frac{1}{2}\mathbb{E}_F[|X - X’|]$ rewards the forecast for spreading its probability mass honestly to reflect true uncertainty (sharpness and calibration). If $F$ collapses to a deterministic point forecast $\hat{y}$, the second term vanishes and $\mathrm{CRPS}$ reduces to the Mean Absolute Error $|\hat{y} - y|$.

2. Horizontal Distance: Wasserstein-1 ($W_1$, Earth Mover’s Distance) and Lipschitz Functions

What if we need to compare two full probability distributions $P$ (with CDF $F$) and $Q$ (with CDF $G$)—for example, the distribution of commuter behavior before vs. after a navigation app changes its routing model?

If we think of $P$ and $Q$ as two piles of earth of total mass $1$, the Wasserstein-1 distance $W_1(P, Q)$ (wasserstein_1d in crates/stochology-core/src/layers.rs) is the minimum physical work (mass $\times$ horizontal distance moved) required to reshape pile $P$ into pile $Q$. In one dimension, Vallender’s theorem proves that this equals the area between their quantile functions (or equivalently between their CDFs):

$$W_1(P, Q) = \int_0^1 \left|F^{-1}(u) - G^{-1}(u)\right| du = \int_{-\infty}^{\infty} |F(x) - G(x)|,dx$$

Why is $W_1$ so central to modern stochastic control and performativity (Chapter 2)? Because of Kantorovich-Rubinstein Duality. A function $h: \mathbb{R} \to \mathbb{R}$ is called $L$-Lipschitz if its slope never exceeds $L$:

$$|h(x) - h(x’)| \le L \cdot |x - x’| \quad \forall x, x’$$

Kantorovich-Rubinstein duality states that for any $L$-Lipschitz loss function $h$, the change in expected loss when the world shifts from distribution $P$ to distribution $Q$ is bounded by:

$$\left|\mathbb{E}_{X \sim P}[h(X)] - \mathbb{E}_{X \sim Q}[h(X)]\right| \le L \cdot W_1(P, Q)$$

This inequality is the mathematical engine of Bridge $L_3 \to L_1$ (Performativity) in Chapter 2: when a deployed model shifts its parameter by $|\theta - \theta’|$, and the world’s distribution shifts by $W_1(\mathcal{D}(\theta), \mathcal{D}(\theta’)) \le \varepsilon |\theta - \theta’|$, Kantorovich-Rubinstein duality tells us whether the feedback loop Between Agent and World will converge to stability or explode!

Step 3 — Auditable Laboratory in Rust and Exercises

All four building blocks of this chapter are implemented in pure Rust in crates/stochology-core/src/layers.rs:

• QuantileSummary::from_samples(&samples): Sorts $n$ observations, computes $\mathbb{E}[X]$, $q_{0.05}$, $q_{0.50}$, $q_{0.95}$, and validates order invariants.
• JensenAudit::evaluate(&samples, f, name): Computes $f(\mathbb{E}[X])$, $\mathbb{E}[f(X)]$, and Jensen’s Gap $\Delta_J = \mathbb{E}[f(X)] - f(\mathbb{E}[X])$, alongside the full QuantileSummary of both $X$ and $f(X)$.
• crps_empirical(&forecast_samples, y_observed): Evaluates the exact $\mathcal{O}(n \log n)$ empirical $\mathrm{CRPS}$ integral from sorted ensemble draws.
• wasserstein_1d(&samples_p, &samples_q): Computes the exact 1D Earth Mover’s distance $W_1(P, Q)$ between two empirical distributions.

Run Lab 1 to watch these primitives operate on synthetic and empirical distributions:

cargo run --release --bin labs -- --layers

Exercises

1. Hand Calculation of Order Statistics and Jensen’s Gap: Given the sample $X = {2, 4, 6, 8, 20}$ and the convex cost function $f(x) = x^2$:
   • Compute $\mathbb{E}[X]$, the median $q_{0.50}$, and the variance $\mathrm{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2$.
   • Show that $f(\mathbb{E}[X]) = 64$, $\mathbb{E}[f(X)] = 104$, and verify the exact identity $\Delta_J = \mathbb{E}[X^2] - (\mathbb{E}[X])^2 = \mathrm{Var}(X) = 40$.
2. Cosine Geometry of Correlation: Construct two centered vectors $\tilde{\mathbf{x}} = (-1, 0, 1)$ and $\tilde{\mathbf{y}} = (1, -2, 1)$ in $\mathbb{R}^3$. Show that their dot product $\langle \tilde{\mathbf{x}}, \tilde{\mathbf{y}} \rangle = 0$ ($\rho_{X,Y} = 0$, zero linear correlation) even though $\tilde{y}_i = 3\tilde{x}_i^2 - 2$ is a deterministic quadratic function of $\tilde{x}_i$! What does this teach about the difference between uncorrelated ($\rho = 0$) and independent?
3. Wasserstein-1 vs. CRPS: Suppose Forecast A predicts a deterministic value $x = 10$ (all mass at $10$), Forecast B predicts a uniform distribution on $[8, 12]$, and the true outcome is $y = 10$. Which forecast has a smaller $\mathrm{CRPS}(F, 10)$? Now suppose the true outcome is $y = 12$: compute $\mathrm{CRPS}(F_A, 12)$ vs. $\mathrm{CRPS}(F_B, 12)$ and explain why ensemble spread protects against tail surprises.