Supervised learning as likelihood maximization

Examples

Instead of a family of functions and a hand-picked loss, we specify a family of conditional probability distributions \(P(y|x, \theta)\) and maximize the likelihood.

The likelihood of the parameters given the training data is

\[ L(\theta) = \prod_{i=1}^n P(y_i|x_i, \theta), \qquad \ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log P(y_i|x_i, \theta). \]

Both have their maximum at the same \(\theta\), because \(\log\) is monotonic. We work with the log for two reasons. Products of small numbers underflow, and sums are easier to differentiate.

Linear regression from this point of view

Take a normal distribution whose mean is a linear function of \(x\) and whose variance is constant,

\[ P(y_i|x_i, \beta_0, \beta_1, \sigma) = \frac{1}{\sqrt{2\pi}\,\sigma} \exp\left(-\frac{(y_i - \beta_0 - \beta_1 x_i)^2}{2\sigma^2}\right). \]

Taking the log and dropping the terms that do not depend on \(\beta\), maximizing \(\ell\) is the same as minimizing \(\sum_i (y_i - \beta_0 - \beta_1 x_i)^2\).

So the maximum likelihood solution of this model is the least squares solution.

The same twenty points as on the previous page, and three parameters now, because the normal distribution has a width as well as a mean. On the left the mean of the model with bands at one and two standard deviations; on the right the log-likelihood over the two parameters of the mean, light where it is high, with the current position marked. Changing σ alone changes the log-likelihood without moving the line, and it does not move the contours either: σ sets how tall the ridge is, not where it lies. That is the previous page’s loss surface seen from the other side, and its ridge sits on the same ellipses.

Choosing a normal distribution with constant variance is the same as choosing the squared error.

The maximum likelihood estimate of the variance is

\[ \hat\sigma^2 = \frac1n\sum_{i=1}^n (y_i - \hat\beta_0 - \hat\beta_1 x_i)^2, \]

with \(n\) in the denominator and not \(n-1\). It is biased, and it is a good exercise to check why.

Checking it numerically

import numpy as np
from scipy.optimize import minimize
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(2)
x = rng.uniform(0, 5, 60)
y = 1.3 + 0.8 * x + rng.normal(0, 0.7, 60)


def neg_log_likelihood(params):
    b0, b1, log_sigma = params
    sigma = np.exp(log_sigma)
    resid = y - b0 - b1 * x
    return float(np.sum(0.5 * np.log(2 * np.pi * sigma**2) + resid**2 / (2 * sigma**2)))


opt = minimize(neg_log_likelihood, x0=[0.0, 0.0, 0.0])
lin = LinearRegression().fit(x.reshape(-1, 1), y)

print(f"maximum likelihood: b0 = {opt.x[0]:.4f}, b1 = {opt.x[1]:.4f}, sigma = {np.exp(opt.x[2]):.4f}")
print(f"least squares:      b0 = {lin.intercept_:.4f}, b1 = {lin.coef_[0]:.4f}")
print(f"sigma by hand:      {np.sqrt(np.mean((y - lin.predict(x.reshape(-1,1)))**2)):.4f}")

1

We optimize \(\log\sigma\) instead of \(\sigma\), so that \(\sigma\) stays positive.

maximum likelihood: b0 = 1.3779, b1 = 0.7625, sigma = 0.6765
least squares:      b0 = 1.3779, b1 = 0.7625
sigma by hand:      0.6765

The two agree to four digits. They are the same estimator.

Why this is useful

Once we think in distributions, the loss follows from the response.

response	distribution	loss
continuous	normal	squared error
two classes	Bernoulli	cross-entropy
\(C\) classes	categorical	cross-entropy
counts	Poisson	Poisson deviance

We do not have to invent a loss. We choose a distribution that fits the response, and the loss comes out of the likelihood. This is the subject of week 3.