Maximum likelihood estimation

In statistics, maximum likelihood estimation (MLE) is a method of estimating the parameters of an assumed probability distribution, given some observed data. This is achieved by maximizing a likelihood function so that, under the assumed statistical model, the observed data is most probable. The point in the parameter space that maximizes the likelihood function is called the maximum likelihood estimate.^[1] The logic of maximum likelihood is both intuitive and flexible, and as such the method has become a dominant means of statistical inference.^[2][3][4]

If the likelihood function is differentiable, the derivative test for finding maxima can be applied. In some cases, the first-order conditions of the likelihood function can be solved analytically; for instance, the ordinary least squares estimator for a linear regression model maximizes the likelihood when the random errors are assumed to have normal distributions with the same variance.^[5]

From the perspective of Bayesian inference, MLE is generally equivalent to maximum a posteriori (MAP) estimation with a prior distribution that is uniform in the region of interest. In frequentist inference, MLE is a special case of an extremum estimator, with the objective function being the likelihood.

We model a set of observations as a random sample from an unknown joint probability distribution which is expressed in terms of a set of parameters. The goal of maximum likelihood estimation is to determine the parameters for which the observed data have the highest joint probability. We write the parameters governing the joint distribution as a vector ${\displaystyle \;\theta =\left[\theta _{1},\,\theta _{2},\,\ldots ,\,\theta _{k}\right]^{\mathsf {T}}\;}$ so that this distribution falls within a parametric family ${\displaystyle \;\{f(\cdot \,;\theta )\mid \theta \in \Theta \}\;,}$ where ${\displaystyle \,\Theta \,}$ is called the parameter space, a finite-dimensional subset of Euclidean space. Evaluating the joint density at the observed data sample ${\displaystyle \;\mathbf {y} =(y_{1},y_{2},\ldots ,y_{n})\;}$ gives a real-valued function, $${\displaystyle {\mathcal {L}}_{n}(\theta )={\mathcal {L}}_{n}(\theta ;\mathbf {y} )=f_{n}(\mathbf {y} ;\theta )\;,}$$ which is called the likelihood function. For independent random variables, ${\displaystyle f_{n}(\mathbf {y} ;\theta )}$ will be the product of univariate density functions: $${\displaystyle f_{n}(\mathbf {y} ;\theta )=\prod _{k=1}^{n}\,f_{k}^{\mathsf {univar}}(y_{k};\theta )~.}$$

The goal of maximum likelihood estimation is to find the values of the model parameters that maximize the likelihood function over the parameter space,^[6] that is: $${\displaystyle {\hat {\theta }}={\underset {\theta \in \Theta }{\operatorname {arg\;max} }}\,{\mathcal {L}}_{n}(\theta \,;\mathbf {y} )~.}$$

Intuitively, this selects the parameter values that make the observed data most probable. The specific value ${\displaystyle ~{\hat {\theta }}={\hat {\theta }}_{n}(\mathbf {y} )\in \Theta ~}$ that maximizes the likelihood function ${\displaystyle \,{\mathcal {L}}_{n}\,}$ is called the maximum likelihood estimate. Further, if the function ${\displaystyle \;{\hat {\theta }}_{n}:\mathbb {R} ^{n}\to \Theta \;}$ so defined is measurable, then it is called the maximum likelihood estimator. It is generally a function defined over the sample space, i.e. taking a given sample as its argument. A sufficient but not necessary condition for its existence is for the likelihood function to be continuous over a parameter space ${\displaystyle \,\Theta \,}$ that is compact.^[7] For an open ${\displaystyle \,\Theta \,}$ the likelihood function may increase without ever reaching a supremum value.

In practice, it is often convenient to work with the natural logarithm of the likelihood function, called the log-likelihood: $${\displaystyle \ell (\theta \,;\mathbf {y} )=\ln {\mathcal {L}}_{n}(\theta \,;\mathbf {y} )~.}$$ Since the logarithm is a monotonic function, the maximum of ${\displaystyle \;\ell (\theta \,;\mathbf {y} )\;}$ occurs at the same value of ${\displaystyle \theta }$ as does the maximum of ${\displaystyle \,{\mathcal {L}}_{n}~.}$^[8] If ${\displaystyle \ell (\theta \,;\mathbf {y} )}$ is differentiable in ${\displaystyle \,\Theta \,,}$ sufficient conditions for the occurrence of a maximum (or a minimum) are $${\displaystyle {\frac {\partial \ell }{\partial \theta _{1}}}=0,\quad {\frac {\partial \ell }{\partial \theta _{2}}}=0,\quad \ldots ,\quad {\frac {\partial \ell }{\partial \theta _{k}}}=0~,}$$ known as the likelihood equations. For some models, these equations can be explicitly solved for ${\displaystyle \,{\widehat {\theta \,}}\,,}$ but in general no closed-form solution to the maximization problem is known or available, and an MLE can only be found via numerical optimization. Another problem is that in finite samples, there may exist multiple roots for the likelihood equations.^[9] Whether the identified root ${\displaystyle \,{\widehat {\theta \,}}\,}$ of the likelihood equations is indeed a (local) maximum depends on whether the matrix of second-order partial and cross-partial derivatives, the so-called Hessian matrix

$${\displaystyle \mathbf {H} \left({\widehat {\theta \,}}\right)={\begin{bmatrix}\left.{\frac {\partial ^{2}\ell }{\partial \theta _{1}^{2}}}\right|_{\theta ={\widehat {\theta \,}}}&\left.{\frac {\partial ^{2}\ell }{\partial \theta _{1}\,\partial \theta _{2}}}\right|_{\theta ={\widehat {\theta \,}}}&\dots &\left.{\frac {\partial ^{2}\ell }{\partial \theta _{1}\,\partial \theta _{k}}}\right|_{\theta ={\widehat {\theta \,}}}\\\left.{\frac {\partial ^{2}\ell }{\partial \theta _{2}\,\partial \theta _{1}}}\right|_{\theta ={\widehat {\theta \,}}}&\left.{\frac {\partial ^{2}\ell }{\partial \theta _{2}^{2}}}\right|_{\theta ={\widehat {\theta \,}}}&\dots &\left.{\frac {\partial ^{2}\ell }{\partial \theta _{2}\,\partial \theta _{k}}}\right|_{\theta ={\widehat {\theta \,}}}\\\vdots &\vdots &\ddots &\vdots \\\left.{\frac {\partial ^{2}\ell }{\partial \theta _{k}\,\partial \theta _{1}}}\right|_{\theta ={\widehat {\theta \,}}}&\left.{\frac {\partial ^{2}\ell }{\partial \theta _{k}\,\partial \theta _{2}}}\right|_{\theta ={\widehat {\theta \,}}}&\dots &\left.{\frac {\partial ^{2}\ell }{\partial \theta _{k}^{2}}}\right|_{\theta ={\widehat {\theta \,}}}\end{bmatrix}}~,}$$

is negative semi-definite at ${\displaystyle {\widehat {\theta \,}}}$, as this indicates local concavity. Conveniently, most common probability distributions – in particular the exponential family – are logarithmically concave.^[10][11]

While the domain of the likelihood function—the parameter space—is generally a finite-dimensional subset of Euclidean space, additional restrictions sometimes need to be incorporated into the estimation process. The parameter space can be expressed as $${\displaystyle \Theta =\left\{\theta :\theta \in \mathbb {R} ^{k},\;h(\theta )=0\right\}~,}$$

where ${\displaystyle \;h(\theta )=\left[h_{1}(\theta ),h_{2}(\theta ),\ldots ,h_{r}(\theta )\right]\;}$ is a vector-valued function mapping ${\displaystyle \,\mathbb {R} ^{k}\,}$ into ${\displaystyle \;\mathbb {R} ^{r}~.}$ Estimating the true parameter ${\displaystyle \theta }$ belonging to ${\displaystyle \Theta }$ then, as a practical matter, means to find the maximum of the likelihood function subject to the constraint ${\displaystyle ~h(\theta )=0~.}$

Theoretically, the most natural approach to this constrained optimization problem is the method of substitution, that is "filling out" the restrictions ${\displaystyle \;h_{1},h_{2},\ldots ,h_{r}\;}$ to a set ${\displaystyle \;h_{1},h_{2},\ldots ,h_{r},h_{r+1},\ldots ,h_{k}\;}$ in such a way that ${\displaystyle \;h^{\ast }=\left[h_{1},h_{2},\ldots ,h_{k}\right]\;}$ is a one-to-one function from ${\displaystyle \mathbb {R} ^{k}}$ to itself, and reparameterize the likelihood function by setting ${\displaystyle \;\phi _{i}=h_{i}(\theta _{1},\theta _{2},\ldots ,\theta _{k})~.}$^[12] Because of the equivariance of the maximum likelihood estimator, the properties of the MLE apply to the restricted estimates also.^[13] For instance, in a multivariate normal distribution the covariance matrix ${\displaystyle \,\Sigma \,}$ must be positive-definite; this restriction can be imposed by replacing ${\displaystyle \;\Sigma =\Gamma ^{\mathsf {T}}\Gamma \;,}$ where ${\displaystyle \Gamma }$ is a real upper triangular matrix and ${\displaystyle \Gamma ^{\mathsf {T}}}$ is its transpose.^[14]

In practice, restrictions are usually imposed using the method of Lagrange which, given the constraints as defined above, leads to the restricted likelihood equations $${\displaystyle {\frac {\partial \ell }{\partial \theta }}-{\frac {\partial h(\theta )^{\mathsf {T}}}{\partial \theta }}\lambda =0}$$ and ${\displaystyle h(\theta )=0\;,}$

where ${\displaystyle ~\lambda =\left[\lambda _{1},\lambda _{2},\ldots ,\lambda _{r}\right]^{\mathsf {T}}~}$ is a column-vector of Lagrange multipliers and ${\displaystyle \;{\frac {\partial h(\theta )^{\mathsf {T}}}{\partial \theta }}\;}$ is the k × r Jacobian matrix of partial derivatives.^[12] Naturally, if the constraints are not binding at the maximum, the Lagrange multipliers should be zero.^[15] This in turn allows for a statistical test of the "validity" of the constraint, known as the Lagrange multiplier test.

Nonparametric maximum likelihood estimation can be performed using the empirical likelihood.

A maximum likelihood estimator is an extremum estimator obtained by maximizing, as a function of θ, the objective function ${\displaystyle {\widehat {\ell \,}}(\theta \,;x)}$. If the data are independent and identically distributed, then we have $${\displaystyle {\widehat {\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 )=\operatorname {\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.^[16] 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.

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 {\widehat {\theta \,}}}$ converges in probability to its true value: $${\displaystyle {\widehat {\theta \,}}_{\mathrm {mle} }\ {\xrightarrow {\text{p}}}\ \theta _{0}.}$$

Under slightly stronger conditions, the estimator converges almost surely (or strongly): $${\displaystyle {\widehat {\theta \,}}_{\mathrm {mle} }\ {\xrightarrow {\text{a.s.}}}\ \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.^[17]

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 {\widehat {\ell \,}}(\theta \mid x)}$ is stochastically equicontinuous.

If one wants to demonstrate that the ML estimator ${\displaystyle {\widehat {\theta \,}}}$ converges to θ₀ almost surely, then a stronger condition of uniform convergence almost surely has to be imposed: $${\displaystyle \sup _{\theta \in \Theta }\left\|\;{\widehat {\ell \,}}(\theta \mid x)-\ell (\theta )\;\right\|\ \xrightarrow {\text{a.s.}} \ 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,^[18] $${\displaystyle {\sqrt {n}}\left({\widehat {\theta \,}}_{\mathrm {mle} }-\theta _{0}\right)\ \xrightarrow {d} \ {\mathcal {N}}\left(0,\,I^{-1}\right)}$$ where I is the Fisher information matrix.

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 {\widehat {\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^[19]

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

It maximizes the so-called profile likelihood:

$${\displaystyle {\bar {L}}(\alpha )=\sup _{\theta :\alpha =g(\theta )}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))^{\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))^{\prime }|={\frac {1}{y}}}$ for ${\displaystyle y>0}$.

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,^[18]

$${\displaystyle {\sqrt {n\,}}\,\left({\widehat {\theta \,}}_{\text{mle}}-\theta _{0}\right)\ \ \xrightarrow {d} \ \ {\mathcal {N}}\left(0,\ {\mathcal {I}}^{-1}\right)~,}$$ where ${\displaystyle ~{\mathcal {I}}~}$ is the Fisher information matrix: $${\displaystyle {\mathcal {I}}_{jk}=\operatorname {\mathbb {E} } \,{\biggl [}\;-{\frac {\partial ^{2}\ln f_{\theta _{0}}(X_{t})}{\partial \theta _{j}\,\partial \theta _{k}}}\;{\biggr ]}~.}$$

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

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)^[20]

$${\displaystyle b_{h}\;\equiv \;\operatorname {\mathbb {E} } {\biggl [}\;\left({\widehat {\theta }}_{\mathrm {mle} }-\theta _{0}\right)_{h}\;{\biggr ]}\;=\;{\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}\;=\;\operatorname {\mathbb {E} } \,{\biggl [}\;{\frac {1}{2}}{\frac {\partial ^{3}\ln f_{\theta _{0}}(X_{t})}{\partial \theta _{i}\;\partial \theta _{j}\;\partial \theta _{k}}}+{\frac {\;\partial \ln f_{\theta _{0}}(X_{t})\;}{\partial \theta _{j}}}\,{\frac {\;\partial ^{2}\ln f_{\theta _{0}}(X_{t})\;}{\partial \theta _{i}\,\partial \theta _{k}}}\;{\biggr ]}~.}$$

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 {\widehat {\theta \,}}_{\text{mle}}^{*}={\widehat {\theta \,}}_{\text{mle}}-{\widehat {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.^[21]

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 \operatorname {\mathbb {P} } (\theta \mid x_{1},x_{2},\ldots ,x_{n})={\frac {f(x_{1},x_{2},\ldots ,x_{n}\mid \theta )\operatorname {\mathbb {P} } (\theta )}{\operatorname {\mathbb {P} } (x_{1},x_{2},\ldots ,x_{n})}}}$$

where ${\displaystyle \operatorname {\mathbb {P} } (\theta )}$ is the prior distribution for the parameter θ and where ${\displaystyle \operatorname {\mathbb {P} } (x_{1},x_{2},\ldots ,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},\ldots ,x_{n}\mid \theta )\operatorname {\mathbb {P} } (\theta )}$ with respect to θ. If we further assume that the prior ${\displaystyle \operatorname {\mathbb {P} } (\theta )}$ is a uniform distribution, the Bayesian estimator is obtained by maximizing the likelihood function ${\displaystyle f(x_{1},x_{2},\ldots ,x_{n}\mid \theta )}$. Thus the Bayesian estimator coincides with the maximum likelihood estimator for a uniform prior distribution ${\displaystyle \operatorname {\mathbb {P} } (\theta )}$.

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.^[22]

Thus, the Bayes Decision Rule is stated as

"decide ${\displaystyle \;w_{1}\;}$ if ${\displaystyle ~\operatorname {\mathbb {P} } (w_{1}|x)\;>\;\operatorname {\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}{\operatorname {arg\;max} }}\;\int _{-\infty }^{\infty }\operatorname {\mathbb {P} } ({\text{ error}}\mid x)\operatorname {\mathbb {P} } (x)\,\operatorname {d} x~}$$ where $${\displaystyle \operatorname {\mathbb {P} } ({\text{ error}}\mid x)=\operatorname {\mathbb {P} } (w_{1}\mid x)~}$$ if we decide ${\displaystyle \;w_{2}\;}$ and ${\displaystyle \;\operatorname {\mathbb {P} } ({\text{ error}}\mid x)=\operatorname {\mathbb {P} } (w_{2}\mid x)\;}$ if we decide ${\displaystyle \;w_{1}\;.}$

By applying Bayes' theorem $${\displaystyle \operatorname {\mathbb {P} } (w_{i}\mid x)={\frac {\operatorname {\mathbb {P} } (x\mid w_{i})\operatorname {\mathbb {P} } (w_{i})}{\operatorname {\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}{\operatorname {arg\;max} }}\,{\bigl [}\,\operatorname {\mathbb {P} } (x\mid w)\,\operatorname {\mathbb {P} } (w)\,{\bigr ]}\;,}$$ where ${\displaystyle h_{\text{Bayes}}}$ is the prediction and ${\displaystyle \;\operatorname {\mathbb {P} } (w)\;}$ is the prior probability.

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}}}$).^[23] 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}}}$.^[24]

Maximum likelihood estimates of parameters can be substituted into expressions for the probability density function,