Maximum likelihood estimation

Properties

A maximum likelihood estimator is an extremum estimator obtained by maximizing, as a function of θ, the objective function ${\displaystyle \hat{\ell }(\theta ;x)}$. If the data are independent and identically distributed, then we have $${\displaystyle \hat{\ell }(\theta ;x)=\sum _{i=1}^n\ln f(x_i\mid \theta ),}$$ this being the sample analogue of the expected log-likelihood ${\displaystyle \ell (\theta )=\mathbb{E}[\ln f(x_i\mid \theta )]}$, where this expectation is taken with respect to the true density.

Maximum-likelihood estimators have no optimum properties for finite samples, in the sense that (when evaluated on finite samples) other estimators may have greater concentration around the true parameter-value. However, like other estimation methods, maximum likelihood estimation possesses a number of attractive limiting properties: As the sample size increases to infinity, sequences of maximum likelihood estimators have these properties:

• Consistency: the sequence of MLEs converges in probability to the value being estimated.
• Equivariance: If ${\displaystyle \hat{\theta }}$ is the maximum likelihood estimator for ${\displaystyle \theta }$, and if ${\displaystyle g(\theta )}$ is a bijective transform of ${\displaystyle \theta }$, then the maximum likelihood estimator for ${\displaystyle \alpha =g(\theta )}$ is ${\displaystyle \hat{\alpha }=g(\hat{\theta })}$. The equivariance property can be generalized to non-bijective transforms, although it applies in that case on the maximum of an induced likelihood function which is not the true likelihood in general.
• Efficiency, i.e. it achieves the Cramér–Rao lower bound when the sample size tends to infinity. This means that no consistent estimator has lower asymptotic mean squared error than the MLE (or other estimators attaining this bound), which also means that MLE has asymptotic normality.
• Second-order efficiency after correction for bias.

Consistency

Under the conditions outlined below, the maximum likelihood estimator is consistent. The consistency means that if the data were generated by ${\displaystyle f(\cdot ;{\theta }_0)}$ and we have a sufficiently large number of observations n, then it is possible to find the value of θ₀ with arbitrary precision. In mathematical terms this means that as n goes to infinity the estimator ${\displaystyle \hat{\theta }}$ converges in probability to its true value:

$${\displaystyle {\hat{\theta }}_{\mathrm{m}\mathrm{l}\mathrm{e}}\stackrel{\text{p}}{\to }{\theta }_0.}$$

Under slightly stronger conditions, the estimator converges almost surely (or strongly):

$${\displaystyle {\hat{\theta }}_{\mathrm{m}\mathrm{l}\mathrm{e}}\stackrel{\text{a.s.}}{\to }{\theta }_0.}$$

In practical applications, data is never generated by ${\displaystyle f(\cdot ;{\theta }_0)}$. Rather, ${\displaystyle f(\cdot ;{\theta }_0)}$ is a model, often in idealized form, of the process generated by the data. It is a common aphorism in statistics that all models are wrong. Thus, true consistency does not occur in practical applications. Nevertheless, consistency is often considered to be a desirable property for an estimator to have.

To establish consistency, the following conditions are sufficient.

1. Identification of the model: $${\displaystyle \theta \ne {\theta }_0\iff f(\cdot \mid \theta )\ne f(\cdot \mid {\theta }_0).}$$ In other words, different parameter values θ correspond to different distributions within the model. If this condition did not hold, there would be some value θ₁ such that θ₀ and θ₁ generate an identical distribution of the observable data. Then we would not be able to distinguish between these two parameters even with an infinite amount of data—these parameters would have been observationally equivalent.
   The identification condition is absolutely necessary for the ML estimator to be consistent. When this condition holds, the limiting likelihood function ℓ(θ|·) has unique global maximum at θ₀.
2. Compactness: the parameter space Θ of the model is compact.

   The identification condition establishes that the log-likelihood has a unique global maximum. Compactness implies that the likelihood cannot approach the maximum value arbitrarily close at some other point (as demonstrated for example in the picture on the right).

   Compactness is only a sufficient condition and not a necessary condition. Compactness can be replaced by some other conditions, such as:

   • both concavity of the log-likelihood function and compactness of some (nonempty) upper level sets of the log-likelihood function, or
   • existence of a compact neighborhood N of θ₀ such that outside of N the log-likelihood function is less than the maximum by at least some ε > 0.
3. Continuity: the function ln f(x | θ) is continuous in θ for almost all values of x: $${\displaystyle \mathbb{P}[\ln f(x\mid \theta )\in C^0(\mathrm{\Theta })]=1.}$$ The continuity here can be replaced with a slightly weaker condition of upper semi-continuity.
4. Dominance: there exists D(x) integrable with respect to the distribution f(x | θ₀) such that $${\displaystyle |\ln f(x\mid \theta )|<D(x)\text{ for all }\theta \in \mathrm{\Theta }.}$$ By the uniform law of large numbers, the dominance condition together with continuity establish the uniform convergence in probability of the log-likelihood: $${\displaystyle \underset{\theta \in \mathrm{\Theta }}{sup}\left|\hat{\ell }(\theta \mid x)-\ell (\theta )\right|\stackrel{\text{p}}{\to }0.}$$

The dominance condition can be employed in the case of i.i.d. observations. In the non-i.i.d. case, the uniform convergence in probability can be checked by showing that the sequence ${\displaystyle \hat{\ell }(\theta \mid x)}$ is stochastically equicontinuous.

If one wants to demonstrate that the ML estimator ${\displaystyle \hat{\theta }}$ converges to θ₀ almost surely, then a stronger condition of uniform convergence almost surely has to be imposed: $${\displaystyle \underset{\theta \in \mathrm{\Theta }}{sup}\Vert \hat{\ell }(\theta \mid x)-\ell (\theta )\Vert \stackrel{\text{a.s.}}{\to }0.}$$

Additionally, if (as assumed above) the data were generated by ${\displaystyle f(\cdot ;{\theta }_0)}$, then under certain conditions, it can also be shown that the maximum likelihood estimator converges in distribution to a normal distribution. Specifically, $${\displaystyle \sqrt{n}\left({\hat{\theta }}_{\mathrm{m}\mathrm{l}\mathrm{e}}-{\theta }_0\right)\stackrel{d}{\to }\mathcal{N}\left(0,I^{-1}\right)}$$ where I is the Fisher information matrix.

Functional invariance

The maximum likelihood estimator selects the parameter value which gives the observed data the largest possible probability (or probability density, in the continuous case). If the parameter consists of a number of components, then we define their separate maximum likelihood estimators, as the corresponding component of the MLE of the complete parameter. Consistent with this, if ${\displaystyle \hat{\theta }}$ is the MLE for ${\displaystyle \theta }$, and if ${\displaystyle g(\theta )}$ is any transformation of ${\displaystyle \theta }$, then the MLE for ${\displaystyle \alpha =g(\theta )}$ is by definition

$${\displaystyle \hat{\alpha }=g(\hat{\theta }).}$$

It maximizes the so-called profile likelihood:

$${\displaystyle \overline{L}(\alpha )=\underset{\theta :\alpha =g(\theta )}{sup}L(\theta ).}$$

The MLE is also equivariant with respect to certain transformations of the data. If ${\displaystyle y=g(x)}$ where ${\displaystyle g}$ is one to one and does not depend on the parameters to be estimated, then the density functions satisfy

$${\displaystyle f_Y(y)=f_X(g^{-1}(y))|(g^{-1}(y){)}^{\mathrm{\prime }}|}$$

and hence the likelihood functions for ${\displaystyle X}$ and ${\displaystyle Y}$ differ only by a factor that does not depend on the model parameters.

For example, the MLE parameters of the log-normal distribution are the same as those of the normal distribution fitted to the logarithm of the data. In fact, in the log-normal case if ${\displaystyle X\sim \mathcal{N}(0,1)}$, then ${\displaystyle Y=g(X)=e^X}$ follows a log-normal distribution. The density of Y follows with ${\displaystyle f_X}$ standard Normal and ${\displaystyle g^{-1}(y)=\log (y)}$, ${\displaystyle |(g^{-1}(y){)}^{\mathrm{\prime }}|=\frac{1}{y}}$ for ${\displaystyle y>0}$.

Efficiency

As assumed above, if the data were generated by ${\displaystyle f(\cdot ;{\theta }_0),}$ then under certain conditions, it can also be shown that the maximum likelihood estimator converges in distribution to a normal distribution. It is √n -consistent and asymptotically efficient, meaning that it reaches the Cramér–Rao bound. Specifically,

$${\displaystyle \sqrt{n}\left({\hat{\theta }}_{\text{mle}}-{\theta }_0\right)\stackrel{d}{\to }\mathcal{N}\left(0,{\mathcal{I}}^{-1}\right),}$$ where ${\displaystyle \mathcal{I}}$ is the Fisher information matrix: $${\displaystyle {\mathcal{I}}_{jk}=\mathbb{E}[-\frac{{\mathrm{\partial }}^2\ln f_{{\theta }_0}(X_t)}{\mathrm{\partial }{\theta }_j\mathrm{\partial }{\theta }_k}].}$$

In particular, it means that the bias of the maximum likelihood estimator is equal to zero up to the order ⁠1/√n ⁠.

Second-order efficiency after correction for bias

However, when we consider the higher-order terms in the expansion of the distribution of this estimator, it turns out that θ_mle has bias of order 1⁄n. This bias is equal to (componentwise)

$${\displaystyle b_h\equiv \mathbb{E}[{\left({\hat{\theta }}_{\mathrm{m}\mathrm{l}\mathrm{e}}-{\theta }_0\right)}_h]=\frac{1}{n}\sum _{i,j,k=1}^m{\mathcal{I}}^{hi}{\mathcal{I}}^{jk}\left(\frac{1}{2}{K}_{ijk}+J_{j,ik}\right)}$$

where ${\displaystyle {\mathcal{I}}^{jk}}$ (with superscripts) denotes the (j,k)-th component of the inverse Fisher information matrix ${\displaystyle {\mathcal{I}}^{-1}}$, and

$${\displaystyle \frac{1}{2}{K}_{ijk}+J_{j,ik}=\mathbb{E}[\frac{1}{2}\frac{{\mathrm{\partial }}^3\ln f_{{\theta }_0}(X_t)}{\mathrm{\partial }{\theta }_i\mathrm{\partial }{\theta }_j\mathrm{\partial }{\theta }_k}+\frac{\mathrm{\partial }\ln f_{{\theta }_0}(X_t)}{\mathrm{\partial }{\theta }_j}\frac{{\mathrm{\partial }}^2\ln f_{{\theta }_0}(X_t)}{\mathrm{\partial }{\theta }_i\mathrm{\partial }{\theta }_k}].}$$

Using these formulae it is possible to estimate the second-order bias of the maximum likelihood estimator, and correct for that bias by subtracting it: $${\displaystyle {\hat{\theta }}_{\text{mle}}^{\ast }={\hat{\theta }}_{\text{mle}}-\hat{b}.}$$ This estimator is unbiased up to the terms of order ⁠1/ n ⁠, and is called the bias-corrected maximum likelihood estimator.

This bias-corrected estimator is second-order efficient (at least within the curved exponential family), meaning that it has minimal mean squared error among all second-order bias-corrected estimators, up to the terms of the order ⁠1/ n² ⁠ . It is possible to continue this process, that is to derive the third-order bias-correction term, and so on. However, the maximum likelihood estimator is not third-order efficient.

Relation to Bayesian inference

A maximum likelihood estimator coincides with the most probable Bayesian estimator given a uniform prior distribution on the parameters. Indeed, the maximum a posteriori estimate is the parameter θ that maximizes the probability of θ given the data, given by Bayes' theorem:

$${\displaystyle \mathbb{P}(\theta \mid x_1,x_2,\dots ,x_n)=\frac{f(x_1,x_2,\dots ,x_n\mid \theta )\mathbb{P}(\theta )}{\mathbb{P}(x_1,x_2,\dots ,x_n)}}$$

where ${\displaystyle \mathbb{P}(\theta )}$ is the prior distribution for the parameter θ and where ${\displaystyle \mathbb{P}(x_1,x_2,\dots ,x_n)}$ is the probability of the data averaged over all parameters. Since the denominator is independent of θ, the Bayesian estimator is obtained by maximizing ${\displaystyle f(x_1,x_2,\dots ,x_n\mid \theta )\mathbb{P}(\theta )}$ with respect to θ. If we further assume that the prior ${\displaystyle \mathbb{P}(\theta )}$ is a uniform distribution, the Bayesian estimator is obtained by maximizing the likelihood function ${\displaystyle f(x_1,x_2,\dots ,x_n\mid \theta )}$. Thus the Bayesian estimator coincides with the maximum likelihood estimator for a uniform prior distribution ${\displaystyle \mathbb{P}(\theta )}$.

Application of maximum-likelihood estimation in Bayes decision theory

In many practical applications in machine learning, maximum-likelihood estimation is used as the model for parameter estimation.

The Bayesian Decision theory is about designing a classifier that minimizes total expected risk, especially, when the costs (the loss function) associated with different decisions are equal, the classifier is minimizing the error over the whole distribution.

Thus, the Bayes Decision Rule is stated as

"decide ${\displaystyle w_1}$ if ${\displaystyle \mathbb{P}(w_1|x)>\mathbb{P}(w_2|x);}$ otherwise decide ${\displaystyle w_2}$"

where ${\displaystyle w_1,w_2}$ are predictions of different classes. From a perspective of minimizing error, it can also be stated as $${\displaystyle w=\underset{w}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbb{P}(\text{ error}\mid x)\mathbb{P}(x)\mathrm{d}x}$$ where $${\displaystyle \mathbb{P}(\text{ error}\mid x)=\mathbb{P}(w_1\mid x)}$$ if we decide ${\displaystyle w_2}$ and ${\displaystyle \mathbb{P}(\text{ error}\mid x)=\mathbb{P}(w_2\mid x)}$ if we decide ${\displaystyle w_1.}$

By applying Bayes' theorem $${\displaystyle \mathbb{P}(w_i\mid x)=\frac{\mathbb{P}(x\mid w_i)\mathbb{P}(w_i)}{\mathbb{P}(x)},}$$ and if we further assume the zero-or-one loss function, which is a same loss for all errors, the Bayes Decision rule can be reformulated as: $${\displaystyle h_{\text{Bayes}}=\underset{w}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}[\mathbb{P}(x\mid w)\mathbb{P}(w)],}$$ where ${\displaystyle h_{\text{Bayes}}}$ is the prediction and ${\displaystyle \mathbb{P}(w)}$ is the prior probability.

Relation to minimizing Kullback–Leibler divergence and cross entropy

Finding ${\displaystyle \hat{\theta }}$ that maximizes the likelihood is asymptotically equivalent to finding the ${\displaystyle \hat{\theta }}$ that defines a probability distribution (${\displaystyle Q_{\hat{\theta }}}$) that has a minimal distance, in terms of Kullback–Leibler divergence, to the real probability distribution from which our data were generated (i.e., generated by ${\displaystyle P_{{\theta }_0}}$). In an ideal world, P and Q are the same (and the only thing unknown is ${\displaystyle \theta }$ that defines P), but even if they are not and the model we use is misspecified, still the MLE will give us the "closest" distribution (within the restriction of a model Q that depends on ${\displaystyle \hat{\theta }}$) to the real distribution ${\displaystyle P_{{\theta }_0}}$.

Prediction bias

Maximum likelihood estimates of parameters can be substituted into expressions for the probability density function, cumulative distribution function, or quantile function, to generate predictions of probabilities or quantiles of out-of-sample events. This method for predicting probabilities is recommended in statistics text-books and actuarial textbooks, and is widely used in the scientific literature. However, maximum likelihood prediction fails to propagate the uncertainty around the maximum likelihood parameter estimates into the prediction. As a result, the predicted probabilities are not well calibrated, and should not be expected to correspond to the frequencies of out-of-sample events. In particular, tail exceedance probabilities and tail exceedance quantiles are typically underestimated, sometimes dramatically. The underestimation is largest when there is little training data, many parameters being estimated, and for the far tail. For cases where this prediction bias is a problem, Bayesian predictions can provide a solution if the prior is chosen so as to reduce or eliminate the bias.

Examples

Discrete uniform distribution

Consider a case where n tickets numbered from 1 to n are placed in a box and one is selected at random (see uniform distribution); thus, the sample size is 1. If n is unknown, then the maximum likelihood estimator ${\displaystyle \hat{n}}$ of n is the number m on the drawn ticket. (The likelihood is 0 for n < m, 1⁄n for n ≥ m, and this is greatest when n = m. Note that the maximum likelihood estimate of n occurs at the lower extreme of possible values {m, m + 1, ...}, rather than somewhere in the "middle" of the range of possible values, which would result in less bias.) The expected value of the number m on the drawn ticket, and therefore the expected value of ${\displaystyle \hat{n}}$, is (n + 1)/2. As a result, with a sample size of 1, the maximum likelihood estimator for n will systematically underestimate n by (n − 1)/2.

Discrete distribution, finite parameter space

Suppose one wishes to determine just how biased an unfair coin is. Call the probability of tossing a 'head' p. The goal then becomes to determine p.

Suppose the coin is tossed 80 times: i.e. the sample might be something like x₁ = H, x₂ = T, ..., x₈₀ = T, and the count of the number of heads "H" is observed.

The probability of tossing tails is 1 − p (so here p is θ above). Suppose the outcome is 49 heads and 31 tails, and suppose the coin was taken from a box containing three coins: one which gives heads with probability p = 1⁄3, one which gives heads with probability p = 1⁄2 and another which gives heads with probability p = 2⁄3. The coins have lost their labels, so which one it was is unknown. Using maximum likelihood estimation, the coin that has the largest likelihood can be found, given the data that were observed. By using the probability mass function of the binomial distribution with sample size equal to 80, number successes equal to 49 but for different values of p (the "probability of success"), the likelihood function (defined below) takes one of three values:

$${\displaystyle \begin{array}{cc}\mathbb{P}[\mathrm{H}=49\mid p=\frac{1}{3}]& =(\genfrac{}{}{0ex}{}{80}{49})(\frac{1}{3}{)}^{49}(1-\frac{1}{3}{)}^{31}\approx 0.000,\\ \mathbb{P}[\mathrm{H}=49\mid p=\frac{1}{2}]& =(\genfrac{}{}{0ex}{}{80}{49})(\frac{1}{2}{)}^{49}(1-\frac{1}{2}{)}^{31}\approx 0.012,\\ \mathbb{P}[\mathrm{H}=49\mid p=\frac{2}{3}]& =(\genfrac{}{}{0ex}{}{80}{49})(\frac{2}{3}{)}^{49}(1-\frac{2}{3}{)}^{31}\approx 0.054.\end{array}}$$

The likelihood is maximized when p = 2⁄3, and so this is the maximum likelihood estimate for p.

Discrete distribution, continuous parameter space

Now suppose that there was only one coin but its p could have been any value 0 ≤ p ≤ 1 . The likelihood function to be maximised is $${\displaystyle L(p)=f_D(\mathrm{H}=49\mid p)=(\genfrac{}{}{0ex}{}{80}{49})p^{49}(1-p{)}^{31},}$$

and the maximisation is over all possible values 0 ≤ p ≤ 1 .

One way to maximize this function is by differentiating with respect to p and setting to zero:

$${\displaystyle \begin{array}{cc}0& =\frac{\mathrm{\partial }}{\mathrm{\partial }p}\left((\genfrac{}{}{0ex}{}{80}{49})p^{49}(1-p{)}^{31}\right),\\ 0& =49p^{48}(1-p{)}^{31}-31p^{49}(1-p{)}^{30}\\ & =p^{48}(1-p{)}^{30}\left[49(1-p)-31p\right]\\ & =p^{48}(1-p{)}^{30}\left[49-80p\right].\end{array}}$$

This is a product of three terms. The first term is 0 when p = 0. The second is 0 when p = 1. The third is zero when p = 49⁄80. The solution that maximizes the likelihood is clearly p = 49⁄80 (since p = 0 and p = 1 result in a likelihood of 0). Thus the maximum likelihood estimator for p is 49⁄80.

This result is easily generalized by substituting a letter such as s in the place of 49 to represent the observed number of 'successes' of our Bernoulli trials, and a letter such as n in the place of 80 to represent the number of Bernoulli trials. Exactly the same calculation yields s⁄n which is the maximum likelihood estimator for any sequence of n Bernoulli trials resulting in s 'successes'.

Continuous distribution, continuous parameter space

For the normal distribution ${\displaystyle \mathcal{N}(\mu ,{\sigma }^2)}$ which has probability density function

$${\displaystyle f(x\mid \mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right),}$$

the corresponding probability density function for a sample of n independent identically distributed normal random variables (the likelihood) is

$${\displaystyle f(x_1,\dots ,x_n\mid \mu ,{\sigma }^2)=\prod _{i=1}^{n}f(x_i\mid \mu ,{\sigma }^2)={\left(\frac{1}{2\pi {\sigma }^2}\right)}^{n/2}\exp \left(-\frac{\sum _{i=1}^n(x_i-\mu {)}^2}{2{\sigma }^2}\right).}$$

This family of distributions has two parameters: θ = (μ, σ); so we maximize the likelihood, ${\displaystyle \mathcal{L}(\mu ,{\sigma }^2)=f(x_1,\dots ,x_n\mid \mu ,{\sigma }^2)}$, over both parameters simultaneously, or if possible, individually.

Since the logarithm function itself is a continuous strictly increasing function over the range of the likelihood, the values which maximize the likelihood will also maximize its logarithm (the log-likelihood itself is not necessarily strictly increasing). The log-likelihood can be written as follows:

$${\displaystyle \log \left(\mathcal{L}(\mu ,{\sigma }^2)\right)=-\frac{n}{2}\log (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{i=1}^n{\left(x_i-\mu \right)}^2}$$

(Note: the log-likelihood is closely related to information entropy and Fisher information.)

We now compute the derivatives of this log-likelihood as follows.

$${\displaystyle \begin{array}{cc}0& =\frac{\mathrm{\partial }}{\mathrm{\partial }\mu }\log \left(\mathcal{L}(\mu ,{\sigma }^2)\right)=0-\frac{-2n(\overline{x}-\mu )}{2{\sigma }^2}.\end{array}}$$ where ${\displaystyle \overline{x}}$ is the sample mean. This is solved by

$${\displaystyle \hat{\mu }=\overline{x}=\sum _{i=1}^n\frac{x_i}{n}.}$$

This is indeed the maximum of the function, since it is the only turning point in μ and the second derivative is strictly less than zero. Its expected value is equal to the parameter μ of the given distribution,

$${\displaystyle \mathbb{E}[\hat{\mu }]=\mu ,}$$

which means that the maximum likelihood estimator ${\displaystyle \hat{\mu }}$ is unbiased.

Similarly we differentiate the log-likelihood with respect to σ and equate to zero:

$${\displaystyle \begin{array}{cc}0& =\frac{\mathrm{\partial }}{\mathrm{\partial }\sigma }\log (\mathcal{L}(\mu ,{\sigma }^2))=-\frac{n}{\sigma }+\frac{1}{{\sigma }^3}\sum _{i=1}^n(x_i-\mu {)}^2.\end{array}}$$

which is solved by

$${\displaystyle {\hat{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n(x_i-\mu {)}^2.}$$

Inserting the estimate ${\displaystyle \mu =\hat{\mu }}$ we obtain

$${\displaystyle {\hat{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n(x_i-\overline{x}{)}^2=\frac{1}{n}\sum _{i=1}^n{x}_i^2-\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n{x}_i{x}_j.}$$

To calculate its expected value, it is convenient to rewrite the expression in terms of zero-mean random variables (statistical error) ${\displaystyle {\delta }_i\equiv \mu -x_i}$. Expressing the estimate in these variables yields

$${\displaystyle {\hat{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n(\mu -{\delta }_i{)}^2-\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n(\mu -{\delta }_i)(\mu -{\delta }_j).}$$

Simplifying the expression above, utilizing the facts that ${\displaystyle \mathbb{E}[{\delta }_i]=0}$ and ${\displaystyle \mathrm{E}[{\delta }_i^2]={\sigma }^2}$, allows us to obtain

$${\displaystyle \mathbb{E}[{\hat{\sigma }}^2]=\frac{n-1}{n}{\sigma }^2.}$$

This means that the estimator ${\displaystyle {\hat{\sigma }}^2}$ is biased for ${\displaystyle {\sigma }^2}$. It can also be shown that ${\displaystyle \hat{\sigma }}$ is biased for ${\displaystyle \sigma }$, but that both ${\displaystyle {\hat{\sigma }}^2}$ and ${\displaystyle \hat{\sigma }}$ are consistent.

Formally we say that the maximum likelihood estimator for ${\displaystyle \theta =(\mu ,{\sigma }^2)}$ is

$${\displaystyle \hat{\theta }=\left(\hat{\mu },{\hat{\sigma }}^2\right).}$$

In this case the MLEs could be obtained individually. In general this may not be the case, and the MLEs would have to be obtained simultaneously.

The normal log-likelihood at its maximum takes a particularly simple form:

$${\displaystyle \log (\mathcal{L}(\hat{\mu },\hat{\sigma }))=\frac{-n}{2}(\log (2\pi {\hat{\sigma }}^2)+1)}$$

This maximum log-likelihood can be shown to be the same for more general least squares, even for non-linear least squares. This is often used in determining likelihood-based approximate confidence intervals and confidence regions, which are generally more accurate than those using the asymptotic normality discussed above.

Non-independent variables

It may be the case that variables are correlated, or more generally, not independent. Two random variables ${\displaystyle y_1}$ and ${\displaystyle y_2}$ are independent only if their joint probability density function is the product of the individual probability density functions, i.e.

$${\displaystyle f(y_1,y_2)=f(y_1)f(y_2)}$$

Suppose one constructs an order-n Gaussian vector out of random variables ${\displaystyle (y_1,\dots ,y_n)}$, where each variable has means given by ${\displaystyle ({\mu }_1,\dots ,{\mu }_n)}$. Furthermore, let the covariance matrix be denoted by ${\displaystyle \mathit{\Sigma }}$. The joint probability density function of these n random variables then follows a multivariate normal distribution given by:

$${\displaystyle f(y_1,\dots ,y_n)=\frac{1}{(2\pi {)}^{n/2}\sqrt{det(\mathit{\Sigma })}}\exp \left(-\frac{1}{2}\left[y_1-{\mu }_1,\dots ,y_n-{\mu }_n\right]{\mathit{\Sigma }}^{-1}{\left[y_1-{\mu }_1,\dots ,y_n-{\mu }_n\right]}^{\mathrm{T}}\right)}$$

In the bivariate case, the joint probability density function is given by:

$${\displaystyle f(y_1,y_2)=\frac{1}{2\pi {\sigma }_1{\sigma }_2\sqrt{1-{\rho }^2}}\exp \left[-\frac{1}{2(1-{\rho }^2)}\left(\frac{(y_1-{\mu }_1{)}^2}{{\sigma }_1^2}-\frac{2\rho (y_1-{\mu }_1)(y_2-{\mu }_2)}{{\sigma }_1{\sigma }_2}+\frac{(y_2-{\mu }_2{)}^2}{{\sigma }_2^2}\right)\right]}$$

In this and other cases where a joint density function exists, the likelihood function is defined as above, in the section "principles," using this density.

Example

${\displaystyle X_1,X_2,\dots ,X_m}$ are counts in cells / boxes 1 up to m; each box has a different probability (think of the boxes being bigger or smaller) and we fix the number of balls that fall to be ${\displaystyle n}$:${\displaystyle x_1+x_2+\cdots +x_m=n}$. The probability of each box is ${\displaystyle p_i}$, with a constraint: ${\displaystyle p_1+p_2+\cdots +p_m=1}$. This is a case in which the ${\displaystyle X_i}$ s are not independent, the joint probability of a vector ${\displaystyle x_1,x_2,\dots ,x_m}$ is called the multinomial and has the form:

$${\displaystyle f(x_1,x_2,\dots ,x_m\mid p_1,p_2,\dots ,p_m)=\frac{n!}{\prod x_i!}\prod p_i^{x_i}=(\genfrac{}{}{0ex}{}{n}{x_1,x_2,\dots ,x_m})p_1^{x_1}{p}_2^{x_2}\cdots p_m^{x_m}}$$

Each box taken separately against all the other boxes is a binomial and this is an extension thereof.

The log-likelihood of this is:

$${\displaystyle \ell (p_1,p_2,\dots ,p_m)=\log n!-\sum _{i=1}^m\log x_i!+\sum _{i=1}^m{x}_i\log p_i}$$

The constraint has to be taken into account and use the Lagrange multipliers:

$${\displaystyle L(p_1,p_2,\dots ,p_m,\lambda )=\ell (p_1,p_2,\dots ,p_m)+\lambda \left(1-\sum _{i=1}^m{p}_i\right)}$$

By posing all the derivatives to be 0, the most natural estimate is derived

$${\displaystyle {\hat{p}}_i=\frac{x_i}{n}}$$

Maximizing log likelihood, with and without constraints, can be an unsolvable problem in closed form, then we have to use iterative procedures.

Iterative procedures

Except for special cases, the likelihood equations $${\displaystyle \frac{\mathrm{\partial }\ell (\theta ;\mathbf{y})}{\mathrm{\partial }\theta }=0}$$

cannot be solved explicitly for an estimator ${\displaystyle \hat{\theta }=\hat{\theta }(\mathbf{y})}$. Instead, they need to be solved iteratively: starting from an initial guess of ${\displaystyle \theta }$ (say ${\displaystyle {\hat{\theta }}_1}$), one seeks to obtain a convergent sequence ${\displaystyle \left\{{\hat{\theta }}_r\right\}}$. Many methods for this kind of optimization problem are available, but the most commonly used ones are algorithms based on an updating formula of the form $${\displaystyle {\hat{\theta }}_{r+1}={\hat{\theta }}_r+{\eta }_r{\mathbf{d}}_r\left(\hat{\theta }\right)}$$

where the vector ${\displaystyle {\mathbf{d}}_r\left(\hat{\theta }\right)}$ indicates the descent direction of the rth "step," and the scalar ${\displaystyle {\eta }_r}$ captures the "step length," also known as the learning rate.

Gradient descent method

(Note: here it is a maximization problem, so the sign before gradient is flipped)

$${\displaystyle {\eta }_r\in {\mathbb{R}}^{+}}$$ that is small enough for convergence and ${\displaystyle {\mathbf{d}}_r\left(\hat{\theta }\right)=\mathrm{\nabla }\ell \left({\hat{\theta }}_r;\mathbf{y}\right)}$

Gradient descent method requires to calculate the gradient at the r-th iteration, but no need to calculate the inverse of second-order derivative, i.e., the Hessian matrix. Therefore, it is computationally faster than Newton–Raphson method.

Newton–Raphson method

$${\displaystyle {\eta }_r=1}$$ and ${\displaystyle {\mathbf{d}}_r\left(\hat{\theta }\right)=-{\mathbf{H}}_r^{-1}\left(\hat{\theta }\right){\mathbf{s}}_r\left(\hat{\theta }\right)}$

where ${\displaystyle {\mathbf{s}}_r(\hat{\theta })}$ is the score and ${\displaystyle {\mathbf{H}}_r^{-1}\left(\hat{\theta }\right)}$ is the inverse of the Hessian matrix of the log-likelihood function, both evaluated the rth iteration. But because the calculation of the Hessian matrix is computationally costly, numerous alternatives have been proposed. The popular Berndt–Hall–Hall–Hausman algorithm approximates the Hessian with the outer product of the expected gradient, such that

$${\displaystyle {\mathbf{d}}_r\left(\hat{\theta }\right)=-{\left[\frac{1}{n}\sum _{t=1}^n\frac{\mathrm{\partial }\ell (\theta ;\mathbf{y})}{\mathrm{\partial }\theta }{\left(\frac{\mathrm{\partial }\ell (\theta ;\mathbf{y})}{\mathrm{\partial }\theta }\right)}^{\mathsf{T}}\right]}^{-1}{\mathbf{s}}_r\left(\hat{\theta }\right)}$$

Quasi-Newton methods

Other quasi-Newton methods use more elaborate secant updates to give approximation of Hessian matrix.

Davidon–Fletcher–Powell formula

DFP formula finds a solution that is symmetric, positive-definite and closest to the current approximate value of second-order derivative: $${\displaystyle {\mathbf{H}}_{k+1}=\left(I-{\gamma }_k{y}_k{s}_k^{\mathsf{T}}\right){\mathbf{H}}_k\left(I-{\gamma }_k{s}_k{y}_k^{\mathsf{T}}\right)+{\gamma }_k{y}_k{y}_k^{\mathsf{T}},}$$

where

$${\displaystyle y_k=\mathrm{\nabla }\ell (x_k+s_k)-\mathrm{\nabla }\ell (x_k),}$$ $${\displaystyle {\gamma }_k=\frac{1}{y_k^{\mathsf{T}}{s}_k},}$$ $${\displaystyle s_k=x_{k+1}-x_k.}$$

Broyden–Fletcher–Goldfarb–Shanno algorithm

BFGS also gives a solution that is symmetric and positive-definite:

$${\displaystyle B_{k+1}=B_k+\frac{y_k{y}_k^{\mathsf{T}}}{y_k^{\mathsf{T}}{s}_k}-\frac{B_k{s}_k{s}_k^{\mathsf{T}}{B}_k^{\mathsf{T}}}{s_k^{\mathsf{T}}{B}_k{s}_k},}$$

where

$${\displaystyle y_k=\mathrm{\nabla }\ell (x_k+s_k)-\mathrm{\nabla }\ell (x_k),}$$ $${\displaystyle s_k=x_{k+1}-x_k.}$$

BFGS method is not guaranteed to converge unless the function has a quadratic Taylor expansion near an optimum. However, BFGS can have acceptable performance even for non-smooth optimization instances

Fisher's scoring

Another popular method is to replace the Hessian with the Fisher information matrix, ${\displaystyle \mathcal{I}(\theta )=\mathbb{E}\left[{\mathbf{H}}_r\left(\hat{\theta }\right)\right]}$, giving us the Fisher scoring algorithm. This procedure is standard in the estimation of many methods, such as generalized linear models.

Although popular, quasi-Newton methods may converge to a stationary point that is not necessarily a local or global maximum, but rather a local minimum or a saddle point. Therefore, it is important to assess the validity of the obtained solution to the likelihood equations, by verifying that the Hessian, evaluated at the solution, is both negative definite and well-conditioned.

History

Early users of maximum likelihood include Carl Friedrich Gauss, Pierre-Simon Laplace, Thorvald N. Thiele, and Francis Ysidro Edgeworth. It was Ronald Fisher however, between 1912 and 1922, who single-handedly created the modern version of the method.

Maximum-likelihood estimation finally transcended heuristic justification in a proof published by Samuel S. Wilks in 1938, now called Wilks' theorem. The theorem shows that the error in the logarithm of likelihood values for estimates from multiple independent observations is asymptotically χ ²-distributed, which enables convenient determination of a confidence region around any estimate of the parameters. The only difficult part of Wilks' proof depends on the expected value of the Fisher information matrix, which is provided by a theorem proven by Fisher. Wilks continued to improve on the generality of the theorem throughout his life, with his most general proof published in 1962.

Reviews of the development of maximum likelihood estimation have been provided by a number of authors.