3. Foundations II: Sampling, Matrix Algebra, Estimation, and Regression

Learning Objectives. By the end of this chapter, you will be able to derive the Law of Large Numbers (LLN) and the Central Limit Theorem (CLT) ($\sigma / \sqrt{n}$); prove and apply Eve’s Law (the Law of Total Variance) to separate within-regime noise from between-parameter ignorance; perform applied matrix algebra for statistical inference, including Positive-Definite Covariance Matrices ($\Sigma$), Condition Numbers ($\kappa(\Sigma)$), and step-by-step Cholesky Factorization ($\Sigma = L L^\top$); derive Ordinary Least Squares ($\hat{\theta}_{\mathrm{OLS}}$), Generalized Least Squares ($\hat{\theta}_{\mathrm{GLS}}$), and Gaussian Log-Likelihood; and interpret Positive-Definite Kernel Functions $k(x, x’)$ as continuous similarity tables that prepare the ground for Gaussian Processes in Chapter 4.

Step 1 — Intuition and the Real-World Paradox

When Ten Thousand Data Points Give You Zero New Information About the Bias

Suppose an artillery officer wants to calibrate the elevation angle $\theta$ of a cannon by firing $n$ shells at a target $1\text{,}000\text{ m}$ away and measuring where they land.

• In a laboratory wind tunnel where the atmospheric drag equation is $100%$ exact, every shell lands at $y_i = f(x_i, \theta) + \varepsilon_i$, where $\varepsilon_i \sim \mathcal{N}(0, \sigma_n^2)$ is independent gusts of wind.
• If a single shell has wind noise $\sigma_n = 20\text{ m}$, averaging $n = 100$ shells shrinks the standard error to $20 / \sqrt{100} = 2\text{ m}$. Averaging $n = 10\text{,}000$ shells shrinks it to $20 / \sqrt{10\text{,}000} = 0.20\text{ m}$.

Encouraged by this $\frac{1}{\sqrt{n}}$ miracle, the officer fires $10\text{,}000$ shells in the field, fits $\hat{\theta}_{\mathrm{OLS}}$ by least squares, and reports that the cannon’s angle is now known to $\pm 0.001^\circ$! Yet the very next morning, when he aims at a ridge $2\text{,}500\text{ m}$ away, every shell misses by $180\text{ m}$.

What went wrong? The officer’s ballistic formula $f(x, \theta)$ omitted the quadratic Coriolis and altitude air-density drop $\delta(x)$. When he minimized the sum of squared errors $\sum (y_i - f(x_i, \theta))^2$ at $x = 1\text{,}000\text{ m}$, the least-squares formula forced the parameter $\hat{\theta}$ to absorb the missing physics $\delta(1\text{,}000)$! Firing $10\text{,}000$ shells only drove the random wind variance $\sigma_n^2 / n$ to zero; it locked the parameter into a false value with near-zero reported uncertainty.

Before we can understand how Chapter 4 diagnoses and cures this pathology—using the Law of Total Variance, Kennedy and O’Hagan’s Gaussian Process discrepancy $\delta(x) \sim \mathcal{GP}(0, k)$, Cholesky factorization $\Sigma = L L^\top$, and Info-Gap robustness—we must build the foundations of classical sampling, linear algebra, least-squares estimation, and kernel covariance matrices from the ground up.

Step 2 — Foundations and Mathematical Formalization

3.1 From Finite Samples to the $\sqrt{n}$ Law and Eve’s Law of Total Variance

Why Averaging $n$ Independent Measurements Shrinks Noise by $\sqrt{n}$

Let $X_1, X_2, \dots, X_n$ be independent and identically distributed (i.i.d.) random variables with mean $\mu = \mathbb{E}[X_i]$ and finite variance $\sigma^2 = \mathrm{Var}(X_i) < \infty$. Consider the sample mean $\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i$:

1. Unbiasedness (by linearity of expectation): $$\mathbb{E}[\bar{X}_n] = \frac{1}{n}\sum_{i=1}^n \mathbb{E}[X_i] = \frac{n\mu}{n} = \mu$$
2. Variance of the Sample Mean (by independence, $\mathrm{Cov}(X_i, X_j) = 0$ for $i \ne j$): Because scaling a random variable by $1/n$ scales its variance by $(1/n)^2 = 1/n^2$, $$\mathrm{Var}(\bar{X}_n) = \mathrm{Var}!\left(\frac{1}{n}\sum_{i=1}^n X_i\right) = \frac{1}{n^2}\sum_{i=1}^n \mathrm{Var}(X_i) = \frac{n\sigma^2}{n^2} = \frac{\sigma^2}{n}$$

Taking the square root gives the Standard Error of the Mean:

$$\mathrm{SE}(\bar{X}_n) = \frac{\sigma}{\sqrt{n}}$$

This establishes two foundational results of classical statistics:

• Law of Large Numbers (LLN): As $n \to \infty$, $\mathrm{Var}(\bar{X}_n) = \sigma^2 / n \to 0$, so $\bar{X}_n$ converges in probability to $\mu$ (provided $\sigma^2 < \infty$—as we will see in Chapter 6, heavy-tailed Pareto processes with $\alpha \le 2$ violate this!).
• Central Limit Theorem (CLT): When finite-variance independent perturbations are summed and scaled by $\sqrt{n}$, the standardized error $Z_n = \frac{\bar{X}_n - \mu}{\sigma / \sqrt{n}}$ converges in distribution to a Standard Gaussian $\mathcal{N}(0, 1)$ with density $\phi(z) = \frac{1}{\sqrt{2\pi}}e^{-z^2/2}$.

Eve’s Law (The Law of Total Variance)

What if the data-generating process depends on an unknown regime or parameter $\theta$? How do we separate the randomness of individual outcomes given $\theta$ from our uncertainty about which $\theta$ is true?

First, recall the Law of Iterated Expectations (Adam’s Law): the overall expectation of $Y$ is the average of its conditional expectations across all possible values of $\theta$:

$$\mathbb{E}[Y] = \mathbb{E}_\theta!\left[\mathbb{E}[Y \mid \theta]\right]$$

Using Adam’s Law, we can decompose the total variance $\mathrm{Var}(Y) = \mathbb{E}[Y^2] - (\mathbb{E}[Y])^2$ by adding and subtracting $\mathbb{E}[(\mathbb{E}[Y \mid \theta])^2]$:

$$\begin{aligned} \mathrm{Var}(Y) &= \mathbb{E}_\theta!\left[\mathbb{E}[Y^2 \mid \theta] - (\mathbb{E}[Y \mid \theta])^2\right] + \left(\mathbb{E}_\theta!\left[(\mathbb{E}[Y \mid \theta])^2\right] - \left(\mathbb{E}_\theta[\mathbb{E}[Y \mid \theta]]\right)^2\right) \ &= \underbrace{\mathbb{E}_\theta!\left[\mathrm{Var}(Y \mid \theta)\right]}_{\text{Average Within-Regime Variance}} + \underbrace{\mathrm{Var}_\theta!\left(\mathbb{E}[Y \mid \theta]\right)}_{\text{Variance Between Regime Means}} \end{aligned}$$

This identity is universally known as Eve’s Law (or the Law of Total Variance). In Chapter 4, it provides the exact mathematical boundary between Aleatoric Uncertainty ($\mathbb{E}[\mathrm{Var}(Y \mid \theta)]$, the irreducible physical scatter inside a known regime) and Epistemic Uncertainty ($\mathrm{Var}(\mathbb{E}[Y \mid \theta])$, our reducible ignorance about the true regime $\theta$).

3.2 Applied Matrix Algebra: Positive-Definite Covariance Matrices and Cholesky Factorization

When we observe a vector of $n$ measurements $\mathbf{y} = (y_1, y_2, \dots, y_n)^\top \in \mathbb{R}^n$, their joint second-order structure is captured by an $n \times n$ Covariance Matrix $\Sigma$:

$$\Sigma_{ij} = \mathrm{Cov}(Y_i, Y_j) = \mathbb{E}!\left[(Y_i - \mu_i)(Y_j - \mu_j)\right]$$

• The diagonal entries $\Sigma_{ii} = \mathrm{Var}(Y_i) \ge 0$ are the individual variances of each measurement.
• The off-diagonal entries $\Sigma_{ij} = \Sigma_{ji}$ measure how strongly measurement $i$ and measurement $j$ co-move. Because $\mathrm{Cov}(Y_i, Y_j) = \mathrm{Cov}(Y_j, Y_i)$, every covariance matrix is symmetric ($\Sigma^\top = \Sigma$).

Why Every Valid Covariance Matrix Is Positive Semi-Definite

Take any weight vector $\mathbf{w} = (w_1, \dots, w_n)^\top \in \mathbb{R}^n$ and form the scalar linear combination $Z = \mathbf{w}^\top \mathbf{Y} = \sum_{i=1}^n w_i Y_i$. Because the variance of any real random variable $Z$ cannot be negative:

$$0 \le \mathrm{Var}(\mathbf{w}^\top \mathbf{Y}) = \sum_{i=1}^n \sum_{j=1}^n w_i w_j \mathrm{Cov}(Y_i, Y_j) = \mathbf{w}^\top \Sigma \mathbf{w}$$

If $\mathbf{w}^\top \Sigma \mathbf{w} > 0$ for every non-zero vector $\mathbf{w} \ne \mathbf{0}$, $\Sigma$ is called Symmetric Positive-Definite (SPD). All of its eigenvalues $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_n > 0$ are strictly positive, and the ratio $\kappa(\Sigma) = \lambda_{\max} / \lambda_{\min}$ is its condition number (measuring how close $\Sigma$ is to becoming singular).

Step-by-Step Cholesky Decomposition ($\Sigma = L L^\top$)

In statistical computing (and inside crates/stochology-core/src/uncertainty.rs), we frequently need to solve a linear system $\Sigma \boldsymbol{\alpha} = \mathbf{r}$ or evaluate the determinant $|\Sigma|$. Directly computing the matrix inverse $\Sigma^{-1}$ is a classic numerical blunder: it wastes flops and amplifies floating-point rounding errors whenever $\kappa(\Sigma)$ is large.

Instead, every Symmetric Positive-Definite matrix $\Sigma$ admits a unique Cholesky Factorization—the matrix equivalent of taking a square root:

$$\Sigma = L L^\top$$

where $L$ is a lower-triangular matrix ($L_{ij} = 0$ for all $j > i$) with strictly positive diagonal entries $L_{ii} > 0$. By writing out the row-by-column product $\Sigma_{ij} = \sum_{k=1}^j L_{ik} L_{jk}$ entry by entry from top-left to bottom-right, we obtain the exact recursive algorithm used in cholesky_decompose:

1. Diagonal entry ($i = j$): $$L_{ii} = \sqrt{\Sigma_{ii} - \sum_{k=1}^{i-1} L_{ik}^2}$$
2. Below-diagonal entry ($i > j$): $$L_{ij} = \frac{1}{L_{jj}}\left(\Sigma_{ij} - \sum_{k=1}^{j-1} L_{ik} L_{jk}\right)$$

Worked $2 \times 2$ Example of Cholesky Factorization

Suppose two nearby sensors have covariance matrix:

$$\Sigma = \begin{bmatrix} 4 & 2 \ 2 & 5 \end{bmatrix}$$

Let us find $L = \begin{bmatrix} L_{11} & 0 \ L_{21} & L_{22} \end{bmatrix}$ such that $L L^\top = \Sigma$:

1. Top-left diagonal: $L_{11} = \sqrt{\Sigma_{11}} = \sqrt{4} = 2$.
2. Bottom-left off-diagonal: $L_{21} = \Sigma_{21} / L_{11} = 2 / 2 = 1$.
3. Bottom-right diagonal: $L_{22} = \sqrt{\Sigma_{22} - L_{21}^2} = \sqrt{5 - 1^2} = \sqrt{4} = 2$.

Check: $\begin{bmatrix} 2 & 0 \ 1 & 2 \end{bmatrix}\begin{bmatrix} 2 & 1 \ 0 & 2 \end{bmatrix} = \begin{bmatrix} 4 & 2 \ 2 & 5 \end{bmatrix} = \Sigma$.

Once we have $L$, two otherwise difficult operations become trivial:

• Solving $\Sigma \boldsymbol{\alpha} = \mathbf{r}$ without matrix inversion: First solve the lower-triangular system $L \mathbf{z} = \mathbf{r}$ from top to bottom (forward substitution), then solve the upper-triangular system $L^\top \boldsymbol{\alpha} = \mathbf{z}$ from bottom to top (backward substitution).
• Computing the Log-Determinant $\ln|\Sigma|$ without overflow: Since $|\Sigma| = |L L^\top| = |L|^2 = (\prod_{i=1}^n L_{ii})^2$, taking logarithms turns the product into a stable sum: $$\ln|\Sigma| = 2 \sum_{i=1}^n \ln(L_{ii})$$

3.3 From Ordinary Least Squares (OLS) to Generalized Least Squares (GLS) and Likelihood

Ordinary Least Squares (OLS) as Orthogonal Projection

Suppose we model a vector of observations $\mathbf{y} \in \mathbb{R}^n$ as a linear combination of $p$ known feature columns collected in an $n \times p$ design matrix $X$:

$$\mathbf{y} = X\boldsymbol{\theta} + \boldsymbol{\varepsilon}$$

How do we find the parameter vector $\hat{\boldsymbol{\theta}}_{\mathrm{OLS}} \in \mathbb{R}^p$ that minimizes the squared Euclidean norm of the residual vector $\mathbf{r}(\boldsymbol{\theta}) = \mathbf{y} - X\boldsymbol{\theta}$?

$$S(\boldsymbol{\theta}) = |\mathbf{y} - X\boldsymbol{\theta}|_2^2 = (\mathbf{y} - X\boldsymbol{\theta})^\top (\mathbf{y} - X\boldsymbol{\theta}) = \mathbf{y}^\top \mathbf{y} - 2\boldsymbol{\theta}^\top X^\top \mathbf{y} + \boldsymbol{\theta}^\top X^\top X \boldsymbol{\theta}$$

Setting the gradient with respect to $\boldsymbol{\theta}$ equal to zero gives the Normal Equations:

$$\nabla_{\boldsymbol{\theta}} S(\boldsymbol{\theta}) = -2 X^\top (\mathbf{y} - X\hat{\boldsymbol{\theta}}_{\mathrm{OLS}}) = \mathbf{0} \implies X^\top X \hat{\boldsymbol{\theta}}_{\mathrm{OLS}} = X^\top \mathbf{y}$$

Notice the geometric meaning of $X^\top(\mathbf{y} - X\hat{\boldsymbol{\theta}}_{\mathrm{OLS}}) = \mathbf{0}$ (shown in Panel 2 of Figure 3.1): the shortest distance from $\mathbf{y}$ to the subspace spanned by the columns of $X$ occurs when the residual vector $\mathbf{r} = \mathbf{y} - X\hat{\boldsymbol{\theta}}_{\mathrm{OLS}}$ is strictly perpendicular (orthogonal) to every column of $X$! Solving for $\hat{\boldsymbol{\theta}}_{\mathrm{OLS}}$ gives:

$$\hat{\boldsymbol{\theta}}_{\mathrm{OLS}} = (X^\top X)^{-1} X^\top \mathbf{y}$$

Why OLS Fails Under Missing Physics $\delta(x)$

What if the true physical world actually generates $\mathbf{y} = X\boldsymbol{\theta}^\ast + \boldsymbol{\delta} + \boldsymbol{\varepsilon}$, where $\boldsymbol{\delta} = (\delta(x_1), \dots, \delta(x_n))^\top$ is a systematic structural discrepancy omitted by our model? Plugging $\mathbf{y}$ into the OLS formula gives:

$$\mathbb{E}[\hat{\boldsymbol{\theta}}_{\mathrm{OLS}}] = (X^\top X)^{-1} X^\top (X\boldsymbol{\theta}^\ast + \boldsymbol{\delta}) = \boldsymbol{\theta}^\ast + \underbrace{(X^\top X)^{-1} X^\top \boldsymbol{\delta}}_{\text{Ghost Parameter Bias!}}$$

Because OLS forces the residual to be orthogonal to $X$, any part of the missing physics $\boldsymbol{\delta}$ that aligns with $X$ is projected directly into $\hat{\boldsymbol{\theta}}_{\mathrm{OLS}}$, distorting the physical parameter even as $n \to \infty$!

Generalized Least Squares (GLS) and Gaussian Log-Likelihood

When the errors between observations are correlated with covariance matrix $\Sigma = \mathrm{Cov}(\boldsymbol{\varepsilon})$ (rather than independent white noise $\sigma_n^2 I$), ordinary Euclidean distance treats redundant, highly correlated points as if they were independent.

Instead, we weight the residuals by the inverse covariance matrix $\Sigma^{-1}$ (which is equivalent to whitening the data by $L^{-1}$, where $\Sigma = L L^\top$), minimizing the Mahalanobis distance $(\mathbf{y} - X\boldsymbol{\theta})^\top \Sigma^{-1}(\mathbf{y} - X\boldsymbol{\theta})$. Setting the gradient to zero yields the Generalized Least Squares (GLS) estimator:

$$\hat{\boldsymbol{\theta}}_{\mathrm{GLS}} = (X^\top \Sigma^{-1} X)^{-1} X^\top \Sigma^{-1} \mathbf{y}$$

Furthermore, if the residual vector $\mathbf{r} = \mathbf{y} - X\boldsymbol{\theta}$ follows a multivariate Gaussian distribution $\mathcal{N}(\mathbf{0}, \Sigma)$, its probability density is:

$$p(\mathbf{y} \mid \boldsymbol{\theta}, \Sigma) = \frac{1}{(2\pi)^{n/2} |\Sigma|^{1/2}} \exp!\left(-\frac{1}{2} \mathbf{r}^\top \Sigma^{-1} \mathbf{r}\right)$$

Taking the natural logarithm gives the Log-Marginal Likelihood objective that we will maximize in Chapter 4:

$$\ln p(\mathbf{y} \mid \boldsymbol{\theta}, \Sigma) = -\frac{1}{2}\underbrace{\mathbf{r}^\top \Sigma^{-1} \mathbf{r}}_{\text{Data Fit Penalty}} - \frac{1}{2}\underbrace{\ln|\Sigma|}_{\text{Complexity Penalty}} - \frac{n}{2}\ln(2\pi)$$

Notice how the Gaussian log-likelihood automatically balances data fit ($\mathbf{r}^\top \Sigma^{-1} \mathbf{r}$, computed via Cholesky triangular solves) against model complexity ($\ln|\Sigma| = 2\sum \ln L_{ii}$)!

3.4 Kernel Functions $k(x, x’)$ as Continuous Similarity Tables

How do we construct the covariance matrix $\Sigma$ when we suspect that a physical simulator has a smooth, unknown structural error $\delta(x)$?

Suppose we don’t know the formula for $\delta(x)$, but we know a simple physical continuity principle: if two input settings $x$ and $x’$ are very close to each other, the model’s structural error $\delta(x)$ and $\delta(x’)$ at those two settings must be strongly correlated. Conversely, if $x$ and $x’$ are far apart, knowing the model’s error at $x$ tells us very little about its error at $x’$.

We encode this smoothness assumption mathematically using a Squared-Exponential (RBF) Covariance Kernel:

$$k(x, x’) = \mathrm{Cov}(\delta(x), \delta(x’)) = \sigma_\delta^2 \exp!\left(-\frac{(x - x’)^2}{2\ell^2}\right)$$

Look at the two hyperparameters of this kernel:

1. Amplitude $\sigma_\delta^2 = k(x, x)$: When $x = x’$, $\exp(0) = 1$, so $k(x, x) = \sigma_\delta^2 = \mathrm{Var}(\delta(x))$ is the marginal variance (vertical scale) of the structural discrepancy.
2. Length-scale $\ell > 0$: Controls how rapidly similarity decays with horizontal distance $|x - x’|$. When $|x - x’| \ll \ell$, $k(x, x’) \approx \sigma_\delta^2$ (errors move together); when $|x - x’| \gg \ell$, $k(x, x’) \to 0$ (errors become independent).

Adding independent sensor measurement noise $\sigma_n^2$ on the diagonal ($i = j$) produces the exact total covariance matrix used by Kennedy and O’Hagan’s Gaussian Process in Chapter 4:

$$\Sigma_{ij} = K_{ij} + \sigma_n^2 \delta_{ij} = \sigma_\delta^2 \exp!\left(-\frac{(x_i - x_j)^2}{2\ell^2}\right) + \sigma_n^2 \mathbf{1}{i = j}$$

Notice that adding the positive “nugget” $\sigma_n^2 I$ shifts every eigenvalue of $K$ upward by $\sigma_n^2 > 0$, guaranteeing that $\Sigma$ is strictly positive-definite and that its Cholesky factorization $\Sigma = L L^\top$ never divides by zero!

Step 3 — Auditable Laboratory in Rust and Exercises

In crates/stochology-core/src/uncertainty.rs, you can inspect the exact zero-dependency Rust implementations of the linear-algebra and estimation primitives introduced in this chapter:

• cholesky_decompose(matrix, n): Computes the lower-triangular factor $L$ such that $\Sigma = L L^\top$, returning None if any diagonal pivot $\le 10^{-12}$.
• cholesky_solve(l, b, n): Executes forward substitution ($L\mathbf{z} = \mathbf{b}$) followed by backward substitution ($L^\top \boldsymbol{\alpha} = \mathbf{z}$) in $\mathcal{O}(n^2)$ time.
• rbf_kernel(x1, x2, sigma_delta, length_scale): Evaluates $k(x_1, x_2)$ and builds the regularized Gram matrix $\Sigma = K + \sigma_n^2 I$.

Run Lab 2 to see OLS bias vs. Cholesky-based GLS/GP calibration in action:

cargo run --release --bin labs -- --uncertainty

Exercises

1. Verifying Eve’s Law by Hand: A factory operates in two equally likely temperature regimes ($\theta \in {0, 1}$, $P(\theta=0)=P(\theta=1)=0.5$). In regime $\theta=0$, output is $Y \mid (\theta=0) \sim \mathcal{N}(10, 4)$ (mean $10$, variance $4$). In regime $\theta=1$, output is $Y \mid (\theta=1) \sim \mathcal{N}(20, 6)$ (mean $20$, variance $6$).
   • Compute the average within-regime (aleatoric) variance $\mathbb{E}[\mathrm{Var}(Y \mid \theta)]$.
   • Compute the between-regime (epistemic) variance $\mathrm{Var}(\mathbb{E}[Y \mid \theta])$.
   • Sum them to find the total marginal variance $\mathrm{Var}(Y)$.
2. Cholesky Forward and Backward Substitution: Using the $2 \times 2$ Cholesky factor $L = \begin{bmatrix} 2 & 0 \ 1 & 2 \end{bmatrix}$ from Worked Example 3.2 (where $\Sigma = \begin{bmatrix} 4 & 2 \ 2 & 5 \end{bmatrix}$) and right-hand side $\mathbf{r} = \begin{bmatrix} 6 \ 11 \end{bmatrix}$:
   • Solve $L\mathbf{z} = \mathbf{r}$ by forward substitution for $\mathbf{z} = (z_1, z_2)^\top$.
   • Solve $L^\top \boldsymbol{\alpha} = \mathbf{z}$ by backward substitution for $\boldsymbol{\alpha} = (\alpha_1, \alpha_2)^\top$, and check that $\Sigma \boldsymbol{\alpha} = \mathbf{r}$ and $\ln|\Sigma| = 2(\ln 2 + \ln 2) = \ln 16$.
3. Ghost Parameter Bias under OLS: Suppose a 1D linear model assumes $y = \theta x$ (through the origin) on inputs $x = (1, 2, 3)^\top$, but the true physical world obeys $y = 2x - 0.5 x^2$ (so $\theta^\ast = 2$ and $\delta(x) = -0.5 x^2$). Compute $\hat{\theta}_{\mathrm{OLS}} = \frac{\sum x_i y_i}{\sum x_i^2}$ by hand and show that OLS estimates $\hat{\theta}_{\mathrm{OLS}} = 2 - \frac{18}{14} \approx 0.714$, missing the true linear slope $\theta^\ast = 2$ by more than $64%$!