Statistical Learning: Essential Mathematics

Taking TMA4268 at NTNU helped me understand several choices I’d been making by habit when training models. I had been more interested in trying neural networks than in studying statistics, but the course gave me better ways to think about model error and uncertainty.

The math underneath the methods

Statistical learning studies how we learn from samples and how well the resulting models might perform on new data. Ideas such as maximum likelihood, hypothesis testing, and the bias-variance decomposition help make the assumptions behind a method explicit.

One way to describe a supervised learning problem is with an input space \(X\), an output space \(Y\), and a distribution \(P\) over pairs of inputs and outputs. We want a function \(f: X \to Y\) that minimizes the expected risk:

\[R(f) = \int_{X \times Y} L(f(x), y) \, dP(x, y)\]

Here \(L\) measures the cost of a prediction. Linear regression, logistic regression, and neural networks can all be studied within this framework, even though they use different functions and losses.

Why bother with statistics

Fitting a model leaves several questions unanswered. We may want to know how uncertain a coefficient is, whether a feature adds useful information, or how much a result changes across samples. Statistical methods give us ways to investigate these questions, provided we check their assumptions.

This linear regression example calculates coefficient estimates, standard errors, and p-values:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

import numpy as np
from scipy import stats

def fit_linear_model(X: np.ndarray, y: np.ndarray) -> dict:
    X = np.column_stack([np.ones(X.shape[0]), X])
    
    beta = np.linalg.inv(X.T @ X) @ X.T @ y
    
    n = X.shape[0]
    y_pred = X @ beta
    residuals = y - y_pred
    mse = np.sum(residuals**2) / (n - X.shape[1])
    var_beta = mse * np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(var_beta))
    
    t_stats = beta / se
    p_values = 2 * (1 - stats.t.cdf(np.abs(t_stats), n - X.shape[1]))
    
    return {
        'coefficients': beta,
        'std_errors': se,
        'p_values': p_values,
    }

The standard errors describe uncertainty in the coefficient estimates under the regression model’s assumptions. The p-values test whether a coefficient differs from zero within that model; they don’t establish that a feature causes the outcome.

The bias-variance decomposition

The bias-variance decomposition was especially useful for understanding squared prediction error:

\[E[(Y - \hat{f}(X))^2] = \text{Var}(\hat{f}(X)) + [\text{Bias}(\hat{f}(X))]^2 + \sigma^2\]

The three terms describe variation across training samples, systematic error in the average prediction, and noise in the outcome. A more flexible model can reduce bias while becoming more sensitive to the particular sample used to train it. This helps explain why reducing training error doesn’t always improve predictions on new data.

Combining ML with statistical rigor

Cross-validation gives several scores, and looking at their spread tells us more than reporting the mean alone. The example below also calculates a simple t-based interval, whose assumptions need care:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

from sklearn.model_selection import KFold
import numpy as np
from scipy import stats

def statistical_cross_validate(model, X, y, n_splits: int = 5):
    kf = KFold(n_splits=n_splits, shuffle=True)
    scores = []
    
    for train_idx, test_idx in kf.split(X):
        X_train, X_test = X[train_idx], X[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        
        model.fit(X_train, y_train)
        scores.append(model.score(X_test, y_test))
    
    mean_score = np.mean(scores)
    ci = stats.t.interval(
        0.95, len(scores) - 1, loc=mean_score, scale=stats.sem(scores)
    )
    
    return {
        'mean_score': mean_score,
        'confidence_interval': ci,
        'std_score': np.std(scores),
    }

The calculation treats the fold scores as independent, even though their training sets overlap. That can make the interval too narrow, so I wouldn’t use it as a reliable 95% confidence interval or as proof that one model is better. Scikit-learn’s example on statistical model comparison explains the dependence and shows a corrected comparison.

Choosing a baseline

The course made me more willing to start with a linear or generalized linear model and examine its errors before adding complexity. That gives me a baseline I can understand. If a neural network improves on it, I still want to know how stable the improvement is and where the model fails.

I now pay more attention to how a result was measured and what assumptions support it. A good score is useful, but I also want to understand how much confidence I can place in it.