Maximum Likelihood View of Linear Regression

In the preceding modules, we developed linear regression from a purely geometric perspective—finding the hyperplane that minimizes the sum of squared residuals. While this geometric view is elegant and computationally convenient, it leaves fundamental questions unanswered:

• Why squared residuals rather than absolute residuals or some other measure?
• How do we quantify uncertainty in our predictions and parameter estimates?
• What assumptions about the data generation process justify our approach?

To answer these questions, we must transition from geometry to probability. The maximum likelihood perspective recasts linear regression as a problem of statistical inference—estimating the parameters of a probabilistic model that describes how the data was generated. At the heart of this probabilistic view lies a critical assumption: the noise in our observations follows a Gaussian (Normal) distribution.

Let's formalize the probabilistic perspective. Instead of viewing regression as curve fitting, we model the data-generating process—how the observed responses $y$ are produced from the inputs $\mathbf{x}$.

The generative model:

We assume that each observation $y_i$ is generated according to:

$$y_i = \mathbf{x}_i^T \boldsymbol{\beta} + \epsilon_i$$

where:

• $\mathbf{x}_i \in \mathbb{R}^p$ is the feature vector for observation $i$
• $\boldsymbol{\beta} \in \mathbb{R}^p$ is the true (unknown) parameter vector we seek to estimate
• $\epsilon_i$ is a random noise term that accounts for all factors not captured by the linear relationship
• $y_i \in \mathbb{R}$ is the observed response

The critical assumption is about the distribution of the noise term $\epsilon_i$. In the Gaussian noise model, we assume:

$$\epsilon_i \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2)$$

Conditional distribution of $y$:

Since $y_i = \mathbf{x}_i^T \boldsymbol{\beta} + \epsilon_i$ and $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$, we can derive the conditional distribution of $y_i$ given $\mathbf{x}_i$:

$$y_i | \mathbf{x}_i \sim \mathcal{N}(\mathbf{x}_i^T \boldsymbol{\beta}, \sigma^2)$$

Or equivalently:

$$p(y_i | \mathbf{x}_i, \boldsymbol{\beta}, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(y_i - \mathbf{x}_i^T \boldsymbol{\beta})^2}{2\sigma^2}\right)$$

This formulation tells us that:

• The mean of $y_i$ is $\mathbf{x}_i^T \boldsymbol{\beta}$ (our prediction)
• The variance of $y_i$ is $\sigma^2$ (the irreducible noise level)
• The response $y_i$ is probabilistically determined—we can only specify its distribution, not its exact value

The Gaussian assumption might seem arbitrary at first. Why not use a uniform distribution? Or a Laplace distribution? Or something else entirely?

The strongest justification comes from the Central Limit Theorem (CLT). In most regression scenarios, the noise term $\epsilon_i$ represents the cumulative effect of many small, independent factors that we haven't (or can't) model explicitly:

• Measurement error in sensors
• Minor variations in environmental conditions
• Small fluctuations in human behavior
• Rounding errors in data collection
• Unexplained biological variability

When the total error is the sum of many small, independent contributions, the CLT guarantees that this sum will be approximately Gaussian—regardless of the distributions of the individual components.

Practical implications:

The CLT justification is powerful because it doesn't require us to know the true distribution of individual error sources. If our noise is genuinely composed of many small effects, Gaussian is a reasonable assumption even if we don't know the exact mechanism.

However, this justification has limits:

1. Finite samples: The CLT is an asymptotic result. With few error sources, the distribution may not be Gaussian.

2. Independence assumption: If error components are correlated, the CLT may not apply.

3. Heavy tails: If some error sources have infinite variance (e.g., Cauchy-distributed outliers), the sum won't be Gaussian.

4. Dominant sources: If one error source dominates, the distribution will reflect that source, not a Gaussian sum.

Beyond the CLT, there are several other compelling reasons to adopt the Gaussian noise model:

1. Maximum Entropy Principle:

Among all probability distributions with a specified mean and variance, the Gaussian distribution has maximum entropy—it assumes the least additional structure beyond those two constraints.

If we know only that:

• The noise has zero mean: $\mathbb{E}[\epsilon] = 0$
• The noise has variance $\sigma^2$: $\mathbb{E}[\epsilon^2] = \sigma^2$

...and we want to choose a distribution that doesn't assume anything else unnecessarily, the Gaussian is the unique answer. It's the "maximally uninformative" distribution subject to our constraints.

2. Analytical Tractability:

The Gaussian distribution leads to closed-form solutions for maximum likelihood estimation. This wasn't historically accidental—Gauss developed much of least-squares theory precisely because Gaussian errors yield tractable mathematics.

3. Stability Under Linear Transformations:

The Gaussian distribution is closed under linear transformations. If $\epsilon \sim \mathcal{N}(0, \sigma^2)$ and $y = \mathbf{x}^T \boldsymbol{\beta} + \epsilon$, then:

$$y | \mathbf{x} \sim \mathcal{N}(\mathbf{x}^T \boldsymbol{\beta}, \sigma^2)$$

This closure property means the entire distribution of $y$ is easily characterized—a property not shared by most other distributions.

4. Connection to Least Squares:

As we'll see in detail, the Gaussian assumption leads directly to the least squares criterion. This provides a probabilistic foundation for a technique that can otherwise seem ad hoc. We're not minimizing squared errors because it's convenient; we're doing it because it's the correct thing to do if the noise is Gaussian.

Let's consolidate our probabilistic model in matrix notation. Given a dataset with $n$ observations and $p$ features:

Data:

• Design matrix $\mathbf{X} \in \mathbb{R}^{n \times p}$ (includes intercept column if needed)
• Response vector $\mathbf{y} \in \mathbb{R}^n$

Model: $$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}$$

where $\boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_n)$—a multivariate Gaussian with zero mean and covariance matrix $\sigma^2 \mathbf{I}_n$.

Equivalently: $$\mathbf{y} | \mathbf{X} \sim \mathcal{N}(\mathbf{X}\boldsymbol{\beta}, \sigma^2 \mathbf{I}_n)$$

The probability density of the entire response vector is:

$$p(\mathbf{y} | \mathbf{X}, \boldsymbol{\beta}, \sigma^2) = \frac{1}{(2\pi\sigma^2)^{n/2}} \exp\left(-\frac{1}{2\sigma^2}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\right)$$

Parameters to estimate:

Our probabilistic model has parameters:

1. $\boldsymbol{\beta} \in \mathbb{R}^p$ — the regression coefficients
2. $\sigma^2 > 0$ — the noise variance

Both are unknown and must be estimated from data. In the maximum likelihood approach (next pages), we find values that maximize the probability of observing our actual data.

Key insight:

Notice that $(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})$ is exactly the residual sum of squares (RSS):

$$\text{RSS}(\boldsymbol{\beta}) = |\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|2^2 = \sum{i=1}^n (y_i - \mathbf{x}_i^T \boldsymbol{\beta})^2$$

This RSS appears in the exponent of the Gaussian density. This is not a coincidence—it's the deep connection between Gaussian noise and least squares that we'll fully exploit.

The Gaussian linear regression model makes several assumptions. Understanding these precisely is crucial for knowing when to trust your results and when to seek alternative methods.

Let's examine each assumption with precision:

Hierarchy of importance:

Not all assumptions are equally critical. In practice:

1. Linearity is fundamental—if violated, the model is systematically wrong
2. Independence matters greatly for standard errors and inference
3. Homoscedasticity affects efficiency and confidence intervals
4. Normality is primarily important for small samples and hypothesis testing

For large samples, the CLT often makes inference approximately valid even without perfect normality. But linearity and independence remain essential.

What happens when the Gaussian noise assumption is violated? The consequences depend on which aspect fails:

1. Non-normality with correct variance:

If the noise is non-Gaussian but still has zero mean and constant variance, the OLS estimator remains:

• Unbiased: $\mathbb{E}[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta}$
• Consistent: $\hat{\boldsymbol{\beta}} \xrightarrow{p} \boldsymbol{\beta}$ as $n \to \infty$
• BLUE: Best Linear Unbiased Estimator (Gauss-Markov still holds!)

The impact is primarily on inference:

• Standard errors may be incorrect for small samples
• Confidence intervals and hypothesis tests may be invalid
• Maximum likelihood interpretation is lost

For large samples, the CLT rescues us—test statistics become approximately normal regardless of error distribution.

2. Heavy tails and outliers:

Gaussian distributions have thin tails—extreme values are exponentially rare. Real data often contains outliers that violate this property.

Quantitatively, for a standard Gaussian:

• Probability of |z| > 2: 4.6%
• Probability of |z| > 3: 0.3%
• Probability of |z| > 4: 0.006%

If your data shows more extreme values than these frequencies suggest, the Gaussian assumption is questionable.

Squared loss (from Gaussian MLE) is especially sensitive to outliers because the loss grows quadratically with the residual. A single observation 10 standard deviations from the mean contributes 100 times as much to the objective as an observation 1 standard deviation away. This can catastrophically distort the fit.

We now arrive at the crucial insight that bridges this entire module. The Gaussian noise assumption doesn't just provide a probabilistic interpretation—it explains why least squares works.

The key result (preview):

Under the model $y_i | \mathbf{x}_i \sim \mathcal{N}(\mathbf{x}_i^T \boldsymbol{\beta}, \sigma^2)$, the maximum likelihood estimator for $\boldsymbol{\beta}$ is exactly the OLS estimator:

$$\hat{\boldsymbol{\beta}}{\text{MLE}} = \hat{\boldsymbol{\beta}}{\text{OLS}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$$

Why does this happen?

Recall the Gaussian density: $$p(\mathbf{y} | \boldsymbol{\beta}, \sigma^2) = \frac{1}{(2\pi\sigma^2)^{n/2}} \exp\left(-\frac{|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|^2}{2\sigma^2}\right)$$

Maximizing this likelihood is equivalent to maximizing the exponent (after taking log), which means minimizing $|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|^2$—the residual sum of squares.

This is profound: Least squares isn't arbitrary; it's the principled choice when noise is Gaussian.

What MLE adds beyond OLS:

The maximum likelihood framework provides additional benefits:

1. Uncertainty quantification: We can derive the sampling distribution of $\hat{\boldsymbol{\beta}}$
2. Variance estimation: We get a principled estimator for $\sigma^2$
3. Model comparison: Likelihood-based criteria (AIC, BIC) enable principled model selection
4. Confidence intervals: Probabilistic statements about parameter values
5. Hypothesis testing: Formal statistical tests for coefficient significance
6. Bayesian extension: A natural starting point for Bayesian regression

These benefits make MLE the foundation for modern statistical inference in regression.

We've established the probabilistic foundation for linear regression. Let's consolidate the key insights:

What's next:

In the next page, we'll derive the log-likelihood function explicitly and show step-by-step how maximizing likelihood leads to the OLS solution. We'll see the beautiful mathematics that connects probability theory to linear algebra, and understand why the Gaussian assumption makes the math tractable while other distributions would not.