Tutorial on Maximum Likelihood Estimation: A Parametric Density Estimation Method

MLE Tutorial in PDF Format

Motivation

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

Also 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 number of HEADS, "H" is observed.

The probability of tossing TAILS is 1 − p. 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. Clearly the probability mass function for this experiment is binomial distribution with
sample size equal to 80, number of successes equal to 49 but different values of p. We have
the following probability mass functions for each of the above mentioned cases:

$ Pr(H = 49 | p = {1}/{3}) = \binom{80}{49}(1/3)^{49}(1 - 1/3)^31 \approx 0.000 $

$ Pr(H = 49 | p = {1}/{2}) = \binom{80}{49}(1/2)^{49}(1 - 1/2)^31 \approx 0.012 $

$ Pr(H = 49 | p = {2}/{3}) = \binom{80}{49}(2/3)^{49}(1 - 2/3)^31 \approx 0.054 $

Based on the above equations, we can conclude that the coin with p = 2 / 3 was more likely
to be picked up for the observations which we were given to begin with.

Definition

The generic situation is that we observe a n-dimensional random vector X with probability
density (or mass) function f(x / θ). It is assumed that θ is a fixed, unknown constant
belonging to the set $ \Theta \subset \mathbb{R}^{n} $.

For $ x \in \mathbb{R}^{n} $, the likelihood function of θ is defined as 

L(θ / x) = f(x / θ)

x is regarded as fixed, and θ is regarded as the variable for L. The log-likelihood function
is defined as l(θ / x) = log L(θ / x)

The Maximum Likelihood Estimate (or MLE) is the value $ \hat{ \theta } = \hat{\theta(x)} \in \Theta $
maximizing L(θ / x), provided it exists:

$ L(\hat{\theta}/(x)) = \underset{\theta}{argmax}[ L(\theta/x) ] $

What is Likelihood function ?

If the probability of an event X dependent on model parameters p is written as
P(X | p)

then we talk about the likelihood

L(p | X)

that is the likelihood of the parameters given the data.

For most sensible models, we will find that certain data are more probable than other data. The aim of maximum likelihood estimation is to find the parameter value(s) that makes the observed data most likely. This is because the likelihood of the parameters given the data is defined to be equal to the probability of the data given the parameters

If we were in the business of making predictions based on a set of solid assumptions, then we would be interested in probabilities - the probability of certain outcomes occurring or not occurring.

However, in the case of data analysis, we have already observed all the data: once they have been observed they are fixed, there is no 'probabilistic' part to them anymore (the word data comes from the Latin word meaning 'given'). We are much more interested in the likelihood of the model parameters that underly the fixed data.

The following is the relation between the likelihood and the probability spaces:

Probability:

                  Knowing Parameters $ \rightarrow $  Prediction of Outcomes

Likelihood:

                  Observation of Data $ \rightarrow $  Estimation of Parameters

Method

Maximum likelihood (ML) estimates need not exist nor be unique. In this section, we show how to compute ML
estimates when they exist and are unique. For computational convenience, the ML estimate
is obtained by maximizing the log-likelihood function, log L(θ / x). This is because the
two functions log L(θ / x) and L(θ / x) are monotonically related to each other
so the same ML estimate is obtained by maximizing either one. Assume that the log-likelihood
function is differentiable, if θ_M'L'E exists, it must satisfy the following partial
differential equation known as the likelihood equation:

$ \frac{d}{d\theta}\left( log L(\theta/x) \right) = 0 $

at θ$ = $θ_M'L'E. This is because maximum or minimum of a continuously differentiable
function implies that its first derivatives vanishes at such points.

The likelihood equation represents a necessary condition for the existence of an MLE estimate.
An additional condition must also be satisfied to ensure that log L(θ / x) is a maximum and not
minimum, since the first derivative cannot reveal this. To be a maximum, the shape of the log-likelihood
function should be convex in the neighborhood of θ_M'L'E. This can be checked by
calculating the second derivatives of the log-likelihoods and showing whether they are all negative
at θ$ = $θ_M'L'E.

$ \frac{d^{2}}{d\theta^{2}}\left( log L(\theta/x) \right) < 0 $

Properties

Some general properties (advantages and disadvantages) of the Maximum Likelihood Estimate are as follows:

• For large data samples (large N) the likelihood function L approaches a Gaussian distribution
• Maximum Likelihood estimates are usually consistent.
• For large N the estimates converge to the true value of the parameters which are estimated.
• Maximum Likelihood Estimates are usually unbiased.
• For all sample sizes the parameter of interest is calculated correctly.
• Maximum Likelihood Estimate is efficient: (the estimates have the smallest variance).
• Maximum Likelihood Estimate is sufficient: (it uses all the information in the observations).
• The solution from the Maximum Likelihood Estimate is unique.

On the other hand, we must know the correct probability distribution for the problem at hand.

Numerical Examples using Maximum Likelihood Estimates

In the following section, we discuss the applications of MLE procedure in estimating unknown
parameters of various common density distributions.

Estimating prior probability using MLE

Consider a two-class classification problem, with classes (ω₁,ω₂) and let
Prob(ω₁) = p and Prob(ω₂) = 1 − p (here p is the unknown parameter). By using
MLE, we can estimate p as follows:

Let the sample be $ \mathcal{D} = (x_{1}, x_{2}, \ldots, x_{N}) $. Let ω_i'j be the class of the feature vector <span class="texhtml" />xj
(N is the sample size and <span class="texhtml" />N₁ is the number of feature vectors belonging to class ω₁). Also assume that
samples$ x_{1}, x_{2} \ldots, x_{N} $ are independent events.

$ Prob(\mathcal{D}/p) = \prod_{j=1}^{N} Prob(\omega_{ij} / p) $
                        $ = p^{N_1} * (1-p)^{N - N_1} $

Please note that $ Prob(\mathcal{D}/p) $ is infinitely differentiable function of $ p $, so the local maxima lies where its derivative is zero.

$ \frac{d}{dp} \left( p^{N_1} * (1-p)^{N - N_1} \right) = 0 $
$ N_1 * p^{N_1 - 1} * (1 - p)^{N - N_1} - (N - N_1) * (1 - p)^{N - N_1 - 1} * p^{N} = 0 $
$ p^{N_1} * (1-p)^{N - N_1 - 1} = 0 $

Solving the above equation for $ p $, we get the following:

$ p^{N_1} * (1-p)^{N - N_1 - 1} = 0 $
$ N_1 * (1 - p) = (N - N_1) * p $

So p is either 0 or 1 by eq. CHANGE HERE and p is $ (N_1/N) $ by eq. CHANGE HERE. Hence this proves that taking the frequencies for probabilities of the feature vectors is optimum and using MLE we showed that p is maximized. Hence the class probabilities are optimum (the likelihood function is maximized using MLE).

Estimating $ \mu $ of a Gaussian distribution when $ \Sigma $ is known

In this section, we estimate the value of $ \mu $ ($ \mu_{MLE} $), when the covariance matrix ($ \Sigma $) in known, for gaussian distribution in n-dimensional feature space.

$ \rho( \mathcal{D} / \mu) = \prod_{j = 1}^{N} \rho ( x_{j} / \mu ) $
log likelihood function is
$ log \rho( \mathcal{D} / \mu) = \sum_{j = 1}^{N} log (\rho ( x_{j} / \mu ) ) $
$ = \sum_{j=1}^{N} log\left[ \frac{1}{ (2\pi)^{\frac{n}{2}} |\Sigma|^{\frac{1}{2}} } exp^{-\frac{(x_j - \mu)^{T} \Sigma^{-1}(x_j - \mu)}{2}} \right] $
$ = \sum_{j=1}^{N} log \left[ \frac{1}{ (2 \pi)^{\frac{n}{2}} |\Sigma|^{\frac{1}{2}} } \right] - \frac{1}{2} \left[(x_j - \mu)^{T} \Sigma^{-1}(x_j - \mu) \right] $

This function is infinitely differentiable function of unknown parameters ($ \mu $)'s. To find the maxima, we set the derivative of this eq. CHANGE HERE to 0.

$ \begin{bmatrix} \frac{d}{d\mu_1} \\ \frac{d}{d\mu_2} \\ \vdots \\ \frac{d}{d\mu_N} \\ \end{bmatrix}_{n \times 1} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ \end{bmatrix}_{n \times 1} $

Differentiating the log likelihood functions yields the following results:

$ \triangledown_{\vec{\mu}} log (\mathcal{D}/\vec{\mu}) = \frac{-1}{2} \sum_{j=1}^{N} \triangledown_{\vec{\mu}} \left[ (x_j - \mu)^{T} \Sigma^{-1}(x_j - \mu) \right] $
$ =\frac{-1}{2} \sum_{j=1}^{N} \begin{bmatrix} \frac{d}{d\mu_1} (x_j - \mu)^{T} \Sigma^{-1}(x_j - \mu) \\ \frac{d}{d\mu_2} (x_j - \mu)^{T} \Sigma^{-1}(x_j - \mu) \\ \vdots \\ \frac{d}{d\mu_N} (x_j - \mu)^{T} \Sigma^{-1}(x_j - \mu) \\ \end{bmatrix}_{n \times 1} $

But we know that,

$ \frac{d}{d\mu_{i}} (x_j - \mu)^{T} \Sigma^{-1}(x_j - \mu) = \left[ \frac{d}{d\mu} (x_j - \mu)^{T} \right] {\Sigma}^{-1}(x_j - \mu) + (x_j - \mu)^{T} \Sigma^{-1} \left[ \frac{d}{d\mu} (x_j - \mu) \right] $
$ = 2 \left[ \frac{d}{d\mu} (x_j - \mu)^{T} \right] {\Sigma}^{-1}(x_j - \mu) $
$ = -2 e_i^{T} {\Sigma}^{-1}(x_j - \mu) $

where $ e_i^{T} = [0,0,0, \ldots, 1, \ldots] $ and 1 is located at position i in the array

$ \triangledown_{\vec{\mu}} log (\mathcal{D}/\vec{\mu}) = \frac{-1}{2} \sum_{j=1}^{N} \begin{bmatrix} -2 e_{1}^{T} {\Sigma}^{-1} (x_j - \mu) \\ -2 e_{2}^{T} {\Sigma}^{-1} (x_j - \mu) \\ \vdots \\ -2 e_{N}^{T} {\Sigma}^{-1} (x_j - \mu) \\ \end{bmatrix}_{n \times 1} $
$ = \sum_{j=1}^{N} \begin{bmatrix} e_{1}^{T} \\ e_{2}^{T} \\ \vdots \\ e_{N}^{T} \\ \end{bmatrix}_{n \times 1} {\Sigma}^{-1} (x_j - \mu) $

Notice that, $ \begin{bmatrix} e_{1}^{T} \\ e_{2}^{T} \\ \vdots \\ e_{N}^{T} \end{bmatrix}_{n \times 1} $ is the identity matrix.

So the above equation reduces to.

$ \triangledown_{\vec{\mu}} log (\mathcal{D}/\vec{\mu}) = \sum_{j=1}^{N} {\Sigma}^{-1} (x_j - \mu) = {\Sigma}^{-1} \sum_{j=1}^{N} (x_j - \mu) $

By solving the equation. CHANGE HERE for $ \mu $.

$ {\Sigma}^{-1} \sum_{j=1}^{N} (x_j - \mu) = 0 $
$ \sum {\Sigma}^{-1} \sum_{j=1}^{N} (x_j - \mu) = 0 $
$ \sum_{j=1}^{N} (x_j - \mu) = 0 $
$ \sum_{j=1}^{N} (x_j) - \sum_{j=1}^{N} (\mu) = 0 $
$ \mu = \frac{1}{N} (\sum_{j=1}^{N}(x_j)) = \mu_{MLE} $

Hence, we proved that using MLE the sample mean is the maximum likelihood estimate of any given sample.

Estimating $ \mu $ and $ \sigma^2 $ for 1-D gaussian distribution using MLE}

In this subsection, we estimate the $ \mu $ and $ \sigma^2 $ for one-dimensional gaussian data. Here $ \theta $ = $ (\theta_1, \theta_2) $ are $ (\mu, \sigma^2) $, which are unknown parameters. We estimate these parameter using the procedure discussed in the section CHANGE HERE.

The log-likelihood function for this case is given by the following equation:
$ log \rho(x_k / \mu, \sigma^2) = log \left[ \frac{1}{\surd{2\pi}\sigma^2}exp(-\frac{x_k - \mu}{2\sigma^2}) \right] $
$ = -\frac{1}{2} log (2\pi\sigma^2) - \frac{1}{2\sigma^2}(x_k - \mu)^2 log \rho(\mathcal{D}/\mu,\sigma^2) $

$ = \prod_{k=1}^{N} log \rho (x_k / \mu,\sigma^2 ) $

Since $ (x_1, x_2, \ldots, x_N) $ are I.I.D's, the density function can be written in product form as follows: $ \rho(\mathcal{D}/\mu,\sigma^2) = \rho (( x_1, x_2, \ldots, x_N ) / \mu,\sigma^2) = \rho(x_1/\mu,\sigma^2)\rho(x_2/\mu,\sigma^2) \ldots \rho(x_N/\mu,\sigma^2) $
$ log \rho(\mathcal{D}/\mu,\sigma^2) = log \prod_{k=1}^{N} \rho (x_k / \mu,\sigma^2 ) $
$ = \sum_{k=1}^{N} \left[ -\frac{1}{2} log (1\pi\sigma^2) - \frac{1}{2\sigma^2} (x_k - \mu)^2 \right] $

Differentiating and setting the eq. ~\ref{eq-1d-mu} to 0, yields the following equations:

$ \triangledown_{\mu,\sigma^2} \rho(\mathcal{D}/\mu,\sigma^2) = \begin{bmatrix} \frac{d}{d\mu} log \rho (\mathcal{D}/\mu,\sigma^2) \\ \frac{d}{d\sigma^2} log \rho (\mathcal{D}/\mu,\sigma^2) \\ \end{bmatrix}_{2 \times 1} $
$ = \begin{bmatrix} \frac{d}{d\mu} (-\frac{N}{2} log (2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{k=1}^{N}(x_k - \mu)) \\ \frac{d}{d\sigma^2} (-\frac{N}{2} log (2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{k=1}^{N}(x_k - \mu))\\ \end{bmatrix}_{2 \times 1} $
$ = \begin{bmatrix} \frac{1}{\sigma^2} \sum_{k=1}^{N}(x_k - \mu) \\ -\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4} \sum_{k=1}^{N}(x_k - \mu)^2\\ \end{bmatrix}_{2 \times 1} $

Solving the eq. ~\ref{eq-1d-mu2} for $ \mu $ and $ \sigma^2 $ yields the following:
$ \frac{1}{\sigma^2} \left( \sum_{k=1}^{N}(x_k - \mu) \right) = 0 $
$ \sum_{k=1}^{N}(x_k - \mu) = 0 $
$ \hat{\mu} = \mu = \frac{1}{N}\sum_{k=1}^{N}(x_k) -\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4} \sum_{k=1}^{N}(x_k - \mu)^2 = 0 $
$ \sum_{k=1}^{N}(x_k - \mu)^2 = N\sigma^2 $
$ \hat{\sigma^2} = \sigma^2 = \frac{ \sum_{k=1}^{N}(x_k - \hat{\mu})^2 }{N} $
which are the mean and standard deviation of the empirical data.

In general for N-dimensional Gaussian data, when X $ \sim $ N($ \mu, \Sigma $) where X $ \in \mathcal{R}^n $, with $ \vec{\mu} $ and $ \Sigma $ are unknown parameters, then MLE for $ \mu $ and $ \Sigma $ are as follows:
$ \hat{\mu} = \frac{1}{N}\sum_{k=1}^{N}(\vec{x_k}) $
$ \hat{\sigma^2} = \frac{ \sum_{k=1}^{N}(\vec{x_k} - \vec{\mu})(\vec{x_k} - \vec{\mu})^T }{N} $

How far does the $ \hat{\mu} $ deviate from the true $ \mu $ when MLE is used.

As discussed in section ~\ref{mu-sigma}, we know that

$ E[\hat{\mu}] = \frac{1}{N}\sum_{k=1}^{N}(\vec{x_k}) = \mu $

Now we compute the expected value of $ (|\hat{\mu} - \mu|)^2 $ as follows:

$ E[ (|\hat{\mu} - \mu|)^2 ] = E[(\hat{\mu} - \mu)(\hat{\mu} - \mu)] $
$ =E[\hat{\mu}\mu - \mu\hat{\mu} - \hat{\mu}\mu + \mu\mu] $
$ =E[\hat{\mu}\hat{\mu}] - 2 E[\mu\hat{\mu}] + E[\mu\mu] $
Substituting E[$ \hat{\mu} $] = $ \mu $ we have $ =E[\hat{\mu}\hat{\mu}] - \mu\mu $
$ =E(\frac{1}{N} \sum_{k=1}^{N}X_k * \sum_{j=1}^{N}X_j) - \mu\mu $
$ =\frac{1}{N^2}\sum_{j,k=1}^{N}E[X_jX_k] - \mu\mu $

Treating $ X_j $ as random variables in the above equations, we have

$ E[ (|\hat{\mu} - \mu|)^2 ] =\frac{1}{N^2}\left[ \sum_{j,k=1}^{N}E(X_j)E(X_k) + \sum_{j,k=1}^{N}E(X_j * X_k) \right] - \mu\mu $
$ =\frac{1}{N^2}\left[ N(N-1)\mu\mu + N*E[X^2] \right] \mu\mu $
$ E[ (|\hat{\mu} - \mu|)^2 ] = \frac{1}{N}E[X^2] - \frac{1}{N} \mu\mu $

But, we know that $ E[|X - \mu|^2] = E[(X - \mu)(X - \mu)] = \sigma^2 $
$ E[X * X] - \mu^2 = \sigma^2 $
$ E[X * X] = \sigma^2 * \mu^2 $

Therefore, we have the following result: $ E[ (|\hat{\mu} - \mu|)^2 ] = \frac{1}{N}(\sigma^2 * \mu^2) - \frac{1}{N}(\mu\mu) $
$ = \frac{1}{N} \sigma^2 $

So the expected value of $ (|\hat{\mu} - \mu|)^2 $ is proportional to the true standard deviation.