Maximum Likelihood and MAP Estimation

Let data $D = \{x_1, \dots, x_n\}$ be drawn i.i.d. from a distribution with density or mass function $p(x \mid \theta)$, where $\theta \in \Theta$ is an unknown parameter vector. The likelihood function is the joint probability of the data viewed as a function of $\theta$ with the data held fixed:

The maximum likelihood estimator is

Because products of many small probabilities underflow numerically and are awkward to differentiate, one almost always works with the log-likelihood:

Since $\log$ is strictly increasing, maximizing $\ell(\theta)$ gives the same maximizer as maximizing $L(\theta)$. In practice one minimizes the negative log-likelihood (NLL), $-\ell(\theta)$, turning estimation into a standard minimization problem compatible with gradient descent and other optimization machinery.

MAP estimation treats $\theta$ itself as random, with a prior density $p(\theta)$ encoding belief before seeing data. Bayes' rule gives the posterior

since $p(D)$ does not depend on $\theta$. The MAP estimator maximizes this posterior:

Comparing the two objectives is the key insight: MAP is MLE plus an extra additive term $\log p(\theta)$. When $p(\theta)$ is a zero-mean Gaussian, $\log p(\theta)$ is (up to a constant) $-\lambda \lVert \theta \rVert_2^2$ — the ridge penalty. When $p(\theta)$ is a Laplace distribution, $\log p(\theta)$ becomes $-\lambda \lVert \theta \rVert_1$ — the lasso penalty. MAP is thus exactly regularized maximum likelihood, with the regularization strength determined by how concentrated the prior is around its mode.

Step 1 — choose a model family. Decide what distribution plausibly generated each data point: Bernoulli for binary labels, Gaussian for continuous noise, categorical/softmax for multi-class labels, Poisson for counts, and so on. This choice encodes assumptions and determines the eventual loss function.

Step 2 — write the log-likelihood. Expand $\ell(\theta) = \sum_i \log p(x_i \mid \theta)$ using the chosen density.

Step 3 — optimize. For simple exponential-family models this can be done in closed form by setting $\nabla_\theta \ell(\theta) = 0$ and solving. For most realistic models (logistic regression, neural networks, mixture models) there is no closed form, and one uses iterative numerical optimization: gradient ascent on $\ell(\theta)$ (equivalently, gradient descent on the NLL), Newton's method, or specialized algorithms like Expectation-Maximization for latent-variable models.

Step 4 (for MAP) — add the prior term. Augment the objective with $\log p(\theta)$ before differentiating or before running gradient descent; this simply adds a penalty gradient $\nabla_\theta \log p(\theta)$ to every update.

Step 5 — check second-order conditions. Verify (via the Hessian, or via convexity of the model class) that the stationary point found is indeed a maximum of the likelihood/posterior, not a saddle point or minimum — important for non-convex models like neural networks, where only local optima are guaranteed. ``` function FIT_MLE(data D = {x_1, ..., x_n}, model family p(x | theta), init theta_0): theta <- theta_0 repeat until convergence:

# negative log-likelihood and its gradient over the whole dataset g <- 0 # gradient accumulator, same shape as theta loss <- 0 for each x_i in D: loss <- loss - log p(x_i | theta) g <- g - grad_theta( log p(x_i | theta) ) theta <- theta - learning_rate g # gradient descent step on NLL return theta function FIT_MAP(data D, model family p(x | theta), prior p(theta), init theta_0): theta <- theta_0 repeat until convergence: g <- - grad_theta( log p(theta) ) # prior contributes a penalty gradient loss <- - log p(theta) for each x_i in D: loss <- loss - log p(x_i | theta) g <- g - grad_theta( log p(x_i | theta) ) theta <- theta - learning_rate g return theta ``` The two routines differ by exactly one term: the prior's negative log-density and its gradient. Everything else — the per-example likelihood gradient, the update rule, the convergence check — is identical. This is the precise sense in which "MAP equals regularized MLE."

Time complexity. For a dataset of $n$ points and a parameter vector of dimension $d$, evaluating the log-likelihood and its gradient once costs $O(nd)$ for most standard models (linear/logistic regression, Gaussians with diagonal covariance), since each of the $n$ terms contributes an $O(d)$ gradient computation. If the model requires a full covariance matrix (multivariate Gaussian MLE), computing and inverting the $d \times d$ covariance costs $O(d^3)$ in addition to the $O(nd^2)$ needed to form it. Iterative optimization (gradient descent, Newton's method) multiplies the per-step cost by the number of iterations $T$ needed for convergence: total cost $O(Tnd)$ for gradient methods, or $O(T(nd^2 + d^3))$ for Newton-type methods that use curvature. Closed-form MLEs (e.g., Gaussian mean and variance, or ordinary least squares via the normal equations) avoid the iteration factor $T$ entirely but may still pay $O(d^3)$ for a matrix inversion or decomposition. Stochastic gradient variants reduce the per-step cost to $O(d)$ per example at the price of noisier updates and more steps to convergence — a classical statistical-computational tradeoff.

Space complexity. Storing the data is $O(nd)$; storing the parameter vector and its gradient is $O(d)$; methods that use second-order information (the Hessian) require $O(d^2)$ space, which becomes prohibitive for large $d$ (e.g., deep networks), motivating first-order or quasi-Newton methods with limited memory.

Correctness of MLE — consistency. Under standard regularity conditions (identifiability of $\theta$, the true parameter $\theta^*$ lying in the interior of $\Theta$, smoothness of $p(x\mid\theta)$), the law of large numbers gives

The right-hand side, as a function of $\theta$, is maximized uniquely at $\theta = \theta^*$. This follows from Gibbs' inequality: for any two densities $p$ and $q$,

with equality iff $p = q$ almost everywhere — precisely the statement that the Kullback–Leibler divergence $D_{\mathrm{KL}}(p \,\|\, q) \ge 0$. So maximizing the empirical log-likelihood is, in the population limit, equivalent to minimizing $D_{\mathrm{KL}}(p(\cdot\mid\theta^*) \,\|\, p(\cdot\mid\theta))$ over $\theta$, whose unique minimizer is $\theta^*$. Hence $\hat\theta_{\text{MLE}} \to \theta^*$ in probability: MLE is consistent. Under further regularity, $\hat\theta_{\text{MLE}}$ is also asymptotically normal and efficient, achieving the Cramér–Rao lower bound on variance among (asymptotically) unbiased estimators as $n \to \infty$.

Correctness of MAP — Bayesian justification. MAP is exactly correct as the mode of the true posterior $p(\theta \mid D)$ by construction of Bayes' rule; no approximation is introduced beyond dropping the normalizing constant $p(D)$, which does not affect the arg max. As $n \to \infty$, the likelihood term dominates a fixed prior (its curvature grows as $O(n)$ while the prior's stays fixed), so $\hat\theta_{\text{MAP}} \to \hat\theta_{\text{MLE}} \to \theta^*$ — the prior's influence vanishes asymptotically, a fact sometimes called the Bernstein–von Mises phenomenon.

Honest tradeoffs.

• MLE can severely overfit when $n$ is small relative to $d$: e.g., MLE variance estimates are biased low, and MLE for logistic regression with linearly separable data diverges to infinity.
• MAP requires choosing a prior family and a strength; a badly mismatched prior introduces bias that does not vanish for finite $n$, and choosing it via the data itself risks circularity (partially addressed by hierarchical/empirical Bayes methods, outside this scope).
• MAP gives a single point estimate and, unlike full Bayesian posterior inference, does not quantify uncertainty about $\theta$ — it reports the mode, not the shape, of the posterior.
• For non-convex models (neural nets), neither MLE nor MAP is guaranteed to find the global optimum; gradient methods only guarantee convergence to a stationary point.

Setup. You flip a possibly biased coin $n = 10$ times and observe $k = 8$ heads. Model each flip as Bernoulli($\theta$), $\theta \in (0,1)$ being the probability of heads.

MLE derivation.

Differentiate and set to zero:

The second derivative $\ell''(\theta) = -k/\theta^2 - (n-k)/(1-\theta)^2 < 0$ confirms a maximum.

MAP derivation with a Beta prior. Suppose you believe the coin is probably close to fair and encode this as a $\mathrm{Beta}(\alpha, \beta)$ prior on $\theta$, with density $p(\theta) \propto \theta^{\alpha-1}(1-\theta)^{\beta-1}$. Take $\alpha = \beta = 5$ (a fairly strong "near 0.5" belief, equivalent to having pre-observed 4 heads and 4 tails). The posterior is

i.e. another Beta distribution (the Beta prior is conjugate to the Bernoulli likelihood). Its mode is

Interpretation table:

This trace shows exactly how MAP behaves as advertised: with a flat (uniform) prior it collapses to MLE; with an informative prior it shrinks the estimate toward the prior mode, and that shrinkage weakens as the sample size grows — the hallmark of regularization. Note the pseudo-counts interpretation of the Beta prior: $\alpha - 1$ acts like extra observed heads and $\beta-1$ like extra observed tails, giving direct intuition for why $\alpha=\beta=5$ pulls a small, unbalanced sample ($8$ heads out of $10$) noticeably toward $0.5$.

Connection to regularized linear regression (second case). For linear regression $y_i = w^\top x_i + \varepsilon_i$ with $\varepsilon_i \sim \mathcal{N}(0,\sigma^2)$, the Gaussian log-likelihood is

so MLE for $w$ is exactly ordinary least squares: minimizing $\sum_i (y_i - w^\top x_i)^2$. Now place a Gaussian prior $w \sim \mathcal{N}(0, \tau^2 I)$. Then

and the MAP objective becomes

which is precisely ridge regression with regularization strength $\lambda = \sigma^2/\tau^2$. A tighter prior (small $\tau^2$, strong belief that weights are near zero) yields a larger $\lambda$ — heavier shrinkage — exactly matching intuition.

Nearly every standard training loss in modern machine learning is a negative log-likelihood in disguise. Training a neural network classifier with cross-entropy loss is exactly MLE under a categorical/softmax output model. Training a regression network with mean squared error is MLE under a Gaussian noise assumption. Weight decay, ubiquitous in deep learning frameworks, is MAP estimation with a Gaussian prior on the weights — an L2 penalty added to the loss. Speech recognition systems historically relied on MLE-trained hidden Markov models; modern language models are trained by maximizing the likelihood of the next token given context. Spam filters, medical diagnosis models, and recommendation systems built on logistic regression or naive Bayes are direct applications of MLE/MAP to Bernoulli and categorical likelihoods. In computational biology, MLE fits population-genetic models to sequence data; in finance, MLE estimates volatility parameters in stochastic models. Bayesian deep learning and probabilistic programming languages (Stan, PyMC, Pyro) build directly on the MAP-as-regularized-MLE view before extending to full posterior approximation via variational inference or MCMC — topics that assume this article's ideas as their starting point.

Ask the tutor: How does maximum likelihood estimation extend to models with latent variables, like Gaussian mixtures, via the Expectation-Maximization algorithm?