# Understanding maximum likelihood estimation

I am told

the method of maximum likelihood says we should use the model that assigns the greatest probability to the data we have observed; formally, the maximum likelihood estimator is found by solving
$$\hat{\theta}= \arg_{\theta}\max\{p(x|\theta)\}$$
where $$p(x|\theta)$$ is called the likelihood function.

Am I correct reading $$p(x|\theta)$$ in English as "the probability of the data given the parameters"? I am confused because we at first seem to be told that MLE is about the probability of the parameters given the data.

[Update]

I am still confused about whether the likelihood function returns probability or probability density. Because Wikipedia says

The likelihood function (often simply called the likelihood) represents the probability of random variable realizations conditional on particular values of the statistical parameters

I am a programmer. When I write a function in code I expect it to return a value of a known type. I want to understand the type that the likelihood function returns.

If the type can be either probability or probability density, why does Wikipedia not make that clear at the start?

• Reading it as "the probability of the data given the parameters" is perfectly accurate. There's no problem speaking of the probability of data, what becomes problematic is to refer to the "probability of parameters". I think the main source of confusion is that this same exact quantity may be viewed as either a probability or a likelihood depending on whether you consider varying x for constant theta, or varying theta for constant x. It's a semantic difference but also has mathematical implications, for instance probabilities must sum to 1 for all x, but likelihoods need not for all theta. Commented Feb 28, 2023 at 23:02
• As I now understand it, likelihood is different from probability in that it can be more than 1. Commented Feb 28, 2023 at 23:17
• Well, p(x|θ) is necessarily less than 1 for any x,θ because it is (also) a probability. But I think maybe what you mean to say is that unlike probabilities, absolute likelihoods are meaningless, only relative likelihoods are meaningful. Therefore, one often doesn't take care to normalize likelihoods in any way (over all x--since we only care about ONE specific x), and just works with whatever proportional value is easiest to compute. Commented Feb 28, 2023 at 23:28
• Commented Mar 2, 2023 at 0:23
• @biohacker: The likelihood $p(x|\theta)$ can be larger than one, because it is the value of a density function. See stats.stackexchange.com/q/485776/7224 and stats.stackexchange.com/a/4221/7224 Commented Mar 2, 2023 at 14:51

Parameters are not random variables but fixed unknowns (at least in the likelihood approach to inference) calibrating the distribution of the observables/observations. Outside a Bayesian setup, it is thus incorrect to talk of a probability distribution on the parameters.

The MLE is the numerical value of the parameters that makes the actual observations the most likely to have occurred $$p(x|\hat\theta) = \max_\theta p(x|\theta)$$ The quantity $$p(x|θ)$$ in the above thus reads as

1. the probability of observing the realisation $$x$$ of the random variable $$X$$ when the parameter value (indexing its distribution) is equal to $$θ$$ in a discrete setting and
2. the density of $$X$$ at the value $$x$$ when the parameter value is equal to $$θ$$ in a continuous setting.

(I would even avoid the term given since this could be interpreted as a conditional probability or density, which does not make sense if $$\theta$$ is not a random variable, i.e., outside a Bayesian framework.)

To quote the very originator of the notion of likelihood, R.A. Fisher in ‘On the mathematical foundations of theoretical statistics’ (1922):

I suggest that we may speak without confusion of the likelihood of one value of p being thrice the likelihood of another (…) likelihood is not here used loosely as a synonym of probability, but simply to express the relative frequencies with which such values of the hypothetical quantity p would in fact yield the observed sample.

The fact that the likelihood function can be defined in a way that includes contributions that are not commensurate (the density and the probability mass) arises from the way in which the likelihood function is defined up to a constant of proportionality, where this "constant" can change with the observation $$x$$, but not with the parameter $$\theta$$.

which reinforces the point that the likelihood function is not a probability density function (for its argument $$\theta$$) and can take numerical values above $$1$$ (or any other bound).

• So is the vertical line symbol ambiguous? Sometimes meaning "given" and sometimes meaning "when" ? As a programmer I feel like I can't get maths to compile! Commented Jan 5, 2021 at 8:25
• I prefer the notation that uses the semicolon ; rather than | when denoting a parameter. Commented Mar 2, 2023 at 3:27
• The Wikipedia page makes a clear distinction between the discrete and continuous cases, between probability mass function and probability density function. Referring to probability in both cases may be a shortcut. The page can be improved as well. Commented Mar 3, 2023 at 21:31
• Thank you @Xi'an for answering my comment about the Wikipedia page ( Sorry I accidentally deleted the comment) Commented Mar 4, 2023 at 3:13

Edit: after revisiting this answer, I realized that some of the language I used previously was not precise enough and could confuse people. In particular, using the term data-generative process when our model is discriminative. I've attempted to fix such issues.

Two years too late but here's my go at it :)

DISCLAIMER: this follows Bishop's amazing book Pattern Recognition and Machine Learning. Some sentences are directly copied (particularly in the deterministic approach section). I have not used quotations because it could break up the flow of the words, as I often make additional statements in between. I'm not sure if that's not good practice, so please let me know if that is not the way to do it.

## Setup

We observe a real-valued input variable $$x$$ and we wish to use this scalar observation to predict the value of a real-valued target variable $$t$$. To this end, we have gathered a set of $$N$$ realizations $$x_i$$ of $$x$$ together with a corresponding set of realizations $$t_i$$ of $$t$$

$$$$\mathbf{x} = (x_1, x_2,...x_N)^T \; , \; \mathbf{t} = (t_1, t_2,...t_N)^T$$$$

In this synthetic curve fitting example, the data is generated by the function $$\text{sin}(2 \pi x)$$, with random noise included in the target values so we can write

$$$$y = \text{sin}(2 \pi x) + N(0, \beta^{-1}).$$$$

Our goal is to exploit the training set in order to make predictions of the value $$\hat{t}$$ of the target variable for some new value $$\hat{x}$$ of the input variable. This involves implicitly trying to find the underlying function which we know to be $$\text{sin}(2 \pi x)$$.

## Deterministic Approach: Least Squares

### Model Setup

In the deterministic--or least squares--approach, we simply consider our output $$t$$ to be a parameterized function of the input variable $$x$$ and weights $$\mathbf{w}$$. For this example, we will consider a hypothesis set of all polynomial functions of the following form:

$$$$y(x, \mathbf{w}) = w_0 + w_1x + w_2x^2 + \dots + w_Mx^M = \sum_{j=0}^M w_jx^j \tag{1}$$$$

where $$M$$ is the order of the polynomial, $$x^j$$ denotes $$x$$ raised to the power of $$j$$, and the polynomial coefficients $$w_0,...,w_M$$ are collectively denoted by the vector $$\mathbf{w}$$. Although $$y(x, \mathbf{w})$$ is a nonlinear function of $$x$$, it is a linear function of the coefficients $$\mathbf{w}$$. Functions that are linear in the unknown weights have special properties and are called linear models.

### Inference

We want to find $$\mathbf{w}$$ so that $$(1)$$ fits the data. To do so, we find a $$\mathbf{w}$$ that minimizes some notion of "error". The fit to the data is constrained by the family of functions we've decided to consider: polynomial functions of order $$M$$.

Our notion of error can be mathematically expressed as an error function that measures the misfit between the function $$y(x, \mathbf{w})$$ for a given value of $$\mathbf{w}$$, and the training data points. A widely used --and simple-- error function is given by the sum-of-squares of the errors between predictions $$y(x_n, \mathbf{w})$$ and the corresponding target values $$t_n$$, so that we minimize

$$$$E(\mathbf{w}) = \frac{1}{2} \sum_{n=1}^N \left( y(x_n, \mathbf{w})- t_n \right)^2 \tag{2}$$$$

where $$\frac{1}{2}$$ is for convenience later on. A Visual representation of the sum-of-squares error function in $$(2)$$ is shown here:

The sum-of-squares error function in $$(2)$$ is computed by taking one half the sum of the squared distances of each data point from the function $$y(x, \mathbf{w})$$. These displacements are shown in red.

Solving for $$\mathbf{w}$$ in this setting is fairly straightforward. Because the error function $$(2)$$ is a quadratic function of the coefficients $$\mathbf{w}$$, its derivative with respect to $$\mathbf{w}$$ will be linear in the elements of $$\mathbf{w}$$. Therefore, the minimization of the error function has a unique solution, denoted $$\mathbf{w}^{\star}$$ which can be found in closed form. The resulting polynomial function is then $$y(x, \mathbf{w}^{\star})$$.

## Maximum Likelihood Approach

In the Least Squares approach we viewed curve fitting purely as an error minimization problem. Now we take a first step in viewing it from a probabilistic perspective so that we can express uncertainty in our predictions.

### Model Setup

As in the least-squares approach we could consider $$x$$ to be related to $$t$$ through a parameterized function $$y(x, \mathbf{w})$$. However, we may not want to make such a strong statement as saying $$t$$ is exactly equal to $$y(x, \mathbf{w})$$. This could be because we think there is noise in the observations $$\mathbf{t}$$, for example due to measurement error. To introduce such uncertainty, we need to place a distribution over the target variable $$t$$. A sensible distributional assumption is to place a Gaussian distribution over $$t$$, with its mean given by the parameterized function $$y(x, \mathbf{w})$$, and its variance being fixed and unknown. This is visualized here:

Illustration of a Gaussian conditional distribution over $$t$$ given $$x$$, where the mean is given by some function of $$x$$, and the variance is fixed.

A perhaps more intuitive perspective is to think about our model describing a process that produces $$t$$ from $$x$$ and parameters $$\mathbf{w}$$. In the least-squares, approach, this process was a deterministic function $$y(x, \mathbf{w})$$. In this section, we extend the process by making the assumption that each $$t$$ is the result of $$y(x, \mathbf{w})$$ and some additive uncertainty, where that additive uncertainty takes the form of a zero-mean Gaussian distribution with unknown variance. This is to say that we are assuming, to have gotten a particular instance of $$t$$:

▪    we are given an instance of $$x$$ 1.
▪    nature used this instance of $$x$$ to get the output of the parameterized function $$y(x, \mathbf{w})$$.
▪    nature sampled from a zero-mean Gaussian, with fixed and unknown variance, and added it
to the output.

This leads to the following model

$$$$p(t | x, \mathbf{w}, \beta) = y(x, \mathbf{w}) + N(0, \beta^{-1}) = N\left( y(x, \mathbf{w}), \beta^{-1} \right) \tag{4}$$$$

where we've used the scaling property of the Gaussian distribution's mean, and we've defined $$\beta = \frac{1}{\sigma^2}$$. $$(4)$$ is referred to as the Gaussian noise model.

### Inference

As the name suggests, in order to use the training data $$( \mathbf{x}, \mathbf{t})$$ to determine the values of the unknown parameters $$\mathbf{w}$$ and $$\beta^{-1}$$, we will search for a setting of $$\mathbf{w}$$ that maximizes the likelihood of the data $$t$$. In other words, we've defined a generative process, and we want to find the correct setting of the parameters $$\mathbf{w}$$ so that the likelihood of our process having created our observed $$\mathbf{t}$$ from our observed $$\mathbf{x}$$ is maximized.

Assuming the data $$(\mathbf{x}, \mathbf{t})$$ were independently sampled from $$(4)$$, the likelihood function is simply the product of each conditional distribution and is evaluted for a particular setting of $$\mathbf{w}$$:

$$$$p(\mathbf{t}|\mathbf{x}, \mathbf{w}, \beta) = \prod_{n=1}^N N(t_n|y(x_n, \mathbf{w}), \beta^{-1}) \tag{5}$$$$

Each time we choose a setting for $$\mathbf{w}$$ and plug it into our model, we are defining a conditional distribution--in particular the one in $$(4)$$. This conditional distribution may agree with the data we have, or it may not. Examples of agreement and disagreement are shown in Figure 3.

A Gaussian noise model shown for a handful of $$x$$, with two different settings for $$\mathbf{w}$$ and $$\beta$$. On the left is a setting of $$\left( \mathbf{w}, \beta \right)$$ that yields a model that disagrees with our observed data. On the right is a setting of $$\left( \mathbf{w}, \beta \right)$$ that yields a model that agrees much better with our observed data. Maximum likelihood looks for the setting of $$\left( \mathbf{w}, \beta\right)$$ that best agrees with our observed data.

We begin by finding the maximum likelihood estimates for $$\mathbf{w}$$. For this example, this amounts to taking the derivative of the likelihood $$(5)$$, setting it equal to zero, and then solving for $$\mathbf{w}$$. So again, finding the maximum likelihood setting for $$\mathbf{w}$$ can be found in closed form. It is common to instead maximize the log likelihood instead of $$(5)$$ for numerical stability and convenience, This which can be written as

$$$$\text{ln} p(\mathbf{t}|\mathbf{x}, \mathbf{w}, \beta^{-1}) = -\frac{\beta}{2} \sum_{n=1}^N \left(y(x_n, \mathbf{w}) - t_n \right)^2 + \frac{N}{2} \text{ln} \beta - \frac{N}{2} \text{ln}(2 \pi) \tag{6}$$$$

For the purpose of taking the derivative of $$(6)$$ with respect to $$\mathbf{w}$$, we can omit the last two terms as they do not depend on $$\mathbf{w}$$. We can also replace the coefficient $$\frac{\beta}{2}$$ with $$\frac{1}{2}$$ since scaling $$(6)$$ by a constant won't chance the location of the maximum with respect to $$\mathbf{w}$$. Lastly, we can equivalently minimize the negative log likelihood. This leaves us with minimizing the following:

$$$$\frac{1}{2} \sum_{n=1}^N \left(y(x_n, \mathbf{w}) - t_n \right)^2 \tag{7}$$$$

and so we see that the sum-of-squares error function has arisen as a consequence of maximizing the likelihood under the assumption of a Gaussian noise distribution.

Once we've found the maximum likelihood estimate for $$\mathbf{w}$$, which we will denote $$\mathbf{w}_{\text{ML}}$$, we can use it to find the setting for the precision parameter $$\beta$$ of the Gaussian conditional distribution. Maximizing $$(6)$$ with respect to $$\beta$$ gives

$$$$\frac{1}{\beta_{\text{ML}}} = \frac{1}{N} \sum_{n=1}^N \left(y(x_n, \mathbf{w}_{\text{ML}}) - t_n \right)^2$$$$

and so we see that the maximum likelihood procedure yields a variance $$\sigma^2$$ being the average squared deviation between the observed data points and the fitted $$y(x, \mathbf{w}_{\text{ML}})$$.

Footnotes:

1 An important distinction is that we are not modelling how we are given $$x$$. To do so would be to model $$p(x)$$, thereby making our model a generative one instead of a discriminative one. Here, we are only modelling $$p(t\vert x)$$, not $$p(t,x) = p(t\vert x) p(x)$$.

• Thank you for your comprehensive go at explaining MLE. However on revisiting my question I realised it was all about what the likelihood function actually measures. Commented Mar 3, 2023 at 9:34