(New page: <br> <center><font size="4"></font> <font size="4">'''Maximum Likelihood Estimators and Examples''' <br> </font> <font size="2">A [https://www.projectrhea.org/learning/slectures.php slec...)
 
 
(7 intermediate revisions by the same user not shown)
Line 12: Line 12:
 
= Outline of the slecture =
 
= Outline of the slecture =
  
* Introduction
+
* Background and Introduction of ML estimator
* Derivation for Maximum Likelihood Estimates (MLE)
+
* Finding Maximum Likelihood Estimators (MLE)
* Examples
+
 
* Summary
 
* Summary
 
* References
 
* References
  
 +
----
 +
----
 +
 +
== Background and Introduction of ML estimator ==
 +
 +
Once we have decided on a model(Probability Distribution), our next step is often to estimate some information from the observed data. There are generally two parametric frameworks for estimating unknown information from data. We will refer to these two general frameworks as the Frequentist and Baysian approaches. One very widely used Frequentist estimator is known as the Maximum Likelihood estimator.
 +
 +
 +
In the frequentist approach, one treats the unknown quantity as a deterministic, but unknown parameter vector, <math> \theta \in \Omega</math>. So for example, after we observe the random vector <math>Y \in \mathbb{R}^{n}</math>, then our objective is to use <math>Y</math> to estimate the unknown scalar or vector <math>\theta</math>. In order to formulate this problem, we will assume that the vector <math>Y</math> has a probability density function given by <math>p_{\theta}(y)</math> where <math>\theta</math> parameterizes a family of density functions for <math>Y</math>. We may then use this family of distributions to determine a function, <math>T : \mathbb{R}^{n} \rightarrow \Omega</math>, that can be used to compute an estimate of the unknown parameter as
 +
 +
<center><math>\hat{\theta} = T(Y)</math></center>
 +
 +
Notice, that since <math>T(Y)</math> is a function of random vector <math>Y</math>, the estimate, <math>\hat{\theta}</math>, is a random vector. The mean of the estimator, <math>\bar{\theta}</math>, can be computed as
 +
 +
<center><math>\bar{\theta} = E_{\theta}[\hat{\theta}] = \int_{\mathbb{R}^{n}} T(y)p_{\theta}(y)dy</math></center>
 +
 +
The difference between the mean of the estimator and the value of the parameter is known as the bias and is given by
 +
 +
<center><math> bias_{\theta} = \bar{\theta} -\theta</math></center>
 +
 +
 +
Similarly, the variance of the estimator is given by
 +
 +
<center><math>var_{\theta} = E_{\theta}[(\hat{\theta} -\bar{\theta})^2]</math></center>
 +
 +
and it is easily shown that the mean squared error (MSE) of the estimate is then given by
 +
 +
<center><math>MSE_{\theta} = E_{\theta}[(\hat{\theta}-\theta)^2] = var_{\theta} + (bias_{\theta})^2</math></center>
 +
 +
 +
 +
Since the bias, variance, and the MSE of the estimator will depend on the specific value of <math>\theta</math>, it is often unclear precisely how to compare the accuracy of different estimators. Even estimators that seem quite poor may produce small or zero error for certain values of <math>\theta</math>. For example, consider the estimator which is fixed to the value <math>\hat{\theta}=1</math>, independent of the data. This would seem to be a very poor estimator, but it has an MSE of 0 when <math>\theta=1</math>.
 +
 +
 +
An estimator is said to be consistent if for all <math>\theta \in \Omega</math>, the MSE of the estimator goes to zero as the number of independent data samples, n, goes to infinity. If an estimator is not consistent, this means that even with arbitrarily large quantities of data, the estimate will not approach the true value of the parameter. Consistency would seem to be the least we would expect of an estimator, but we will later see that even some very intuitive estimators are not always consistent.
 +
 +
 +
Ideally, it would be best if one could select an estimator which has uniformly low bias and variance for all values of <math>\theta</math>. This is not always possible, but when it is we have names for such estimators. For example, <math>\hat{\theta}</math> is said to be an unbiased estimator if for all values of <math>\theta</math> the bias is zero, i.e. <math>\theta = \bar{\theta}</math>. If in addition, for all values of <math>\theta</math>, the variance of an estimator is less than or equal to that of all other unbiased estimators, then we say that the estimator is uniformly minimum variance unbiased(UMVU) estimator.
 +
 +
 +
There are many excellent estimators that have been proposed through the years for many different types of problems. However, one very widely used Frequentist estimator is known as the maximum likelihood(ML) estimator given by
 +
 +
<center><math>\hat{\theta} = \arg\max_{\theta \in \Omega} p_{\theta}(Y)</math></center>
 +
 +
<center><math> = \arg\max_{\theta \in \Omega} \log p_{\theta}(Y)</math></center>
 +
 +
where the notation "argmax" denotes the value of the argument that achieves the global maximum of the function. Notice that these formulas for the ML estimate actually use the random variable <math>Y</math> as an argument to the density function <math>p_{\theta}(y)</math>. This implies that <math>\hat{\theta}</math> is function of <math>Y</math>, which in turn means that <math>\hat{\theta}</math> is a random variable.
 +
 +
 +
When the density function, <math>p_{\theta}(y)</math>, is a continuously differentiable function of <math>\theta</math>, then a necessary condition when computing the ML estimate is that the gradient of the likelihood is zero.
 +
 +
<center><math>\bigtriangledown_{\theta} p_{\theta}(Y)|_{\theta=\hat{\theta}} = 0</math></center>
 +
 +
 +
While the ML estimate is generally not unbiased, it does have a number of desirable properties. It can be shown that under relatively mild technical conditions, the ML estimate is both consistent and asymptotically efficient. This means that the ML estimate attains the Cramer-Rao bound asymptotically as the number of independent data samples,n, goes to infinity. Since this Cramer-Rao bound is a lower bound on the variance of an unbiased estimator, this means that the ML estimate is about as good as any unbiased estimator can be. So if one has no idea what values of <math>\theta</math> may occur, and needs to guarantee good performance for all cases, then the ML estimator is usually a good choice.
 +
 +
 +
----
 +
----
 +
 +
== Finding the ML estimators ==
 +
 +
=Example 1: 1D Gaussuan, <math>\theta</math> unknown, <math>\sigma^2</math> known=
 +
Let <math>\{Y_k\}_{k=1}^n</math> be i.i.d. random variables with Gaussian distribution <math>N(\theta,1)</math> and unknown mean parameter <math>\theta</math>. For this case, the logarithm of the density function is given by:
 +
 +
<center><math>\log p_{\theta}(Y) = -\frac{1}{2\sigma^2}\sum_{k=1}^{n}(Y_k-\theta)^2 - \frac{n}{2} \log (2\pi\sigma^2) </math></center>
 +
 +
 +
Local max are where gradient is zero. Differentiating the log likelihood results in the following expression.
 +
 +
<center><math>\bigtriangledown_{\theta} \log p_{\theta}(Y) = \frac{d\log p_{\theta}(Y)}{d\theta} |_{\theta = \hat{\theta}} = \frac{1}{\sigma^2}\sum_{k=1}^{n}(Y_k - \hat{\theta}) = 0 </math></center>
 +
 +
From this we obtain the ML estimate for <math>\theta</math>. So
 +
 +
<center><math>\hat{\theta}=\frac{1}{n}\sum_{k=1}^{n} y_k</math></center>
 +
 +
Notice that, by the following argument, this particular ML estimate is unbiased.
 +
 +
<center><math>E_{\theta}[\hat{\theta}] = E_{\theta}[\frac{1}{n}\sum_{k=1}^{n}Y_k] = \frac{1}{n}\sum_{k=1}^{n} E_{\theta}[Y_k] = \frac{1}{n} \sum_{k=1}^{n} \theta = \theta </math></center>
 +
 +
Moreover, for this special case, it can be shown that the ML estimator is also the UMVU estimator.
 +
 +
 +
 +
 +
=Example 2: 1D Gaussuan, <math>(\mu,\sigma^2)</math> unknown=
 +
 +
Let <math>\{Y_k\}_{k=1}^n</math> be i.i.d. random variables with Gaussian distribution <math>N(\mu,\sigma^2)</math> and unknown mean parameter <math>\theta</math>. For this case, the density function is giving by:
 +
 +
<center><math>p_{\theta}(Y) = (2\pi\sigma^2)^{-\frac{n}{2}} \exp(- \frac{1}{2\sigma^2}\sum_{k=1}^{n}(Y_k - \mu)) </math></center>
 +
 +
The logarithm of the density function is given by:
 +
 +
<center><math>\log p_{\theta}(Y) = -\frac{1}{2\sigma^2}\sum_{k=1}^{n}(Y_k-\mu)^2 - \frac{n}{2} \log (2\pi\sigma^2) </math></center>
 +
 +
Similarly, differentiating the log likelihood results in the following expressions.
 +
 +
<math> \bigtriangledown_{\mu} \log p_{\theta}(Y) =  \bigtriangledown_{\mu} (-\frac{n}{2}\log 2\pi-\frac{n}{2}\log \sigma^2 - \frac{1}{2\sigma^2} \sum_{k=1}^{n} (Y_k-\mu)^2)</math>
 +
 +
<math>  = -\frac{1}{2\sigma^2}\cdot 2\cdot (-1) \sum_{k=1}^{n}(Y_k-\hat{\mu}) </math>
 +
 +
<math> = \frac{1}{\sigma^2} \sum_{k=1}^{n}(Y_k-\hat{\mu})</math>
 +
 +
<math> = \frac{1}{\sigma^2}(\sum_{k=1}^{n}Y_k  -\hat{\mu}n) = 0  --- 1</math>
 +
 +
 +
 +
<math> \bigtriangledown_{\sigma^2} \log p_{\theta}(Y) =  \bigtriangledown_{\sigma^2} (-\frac{n}{2}\log 2\pi-\frac{n}{2}\log \sigma^2 - \frac{1}{2\sigma^2} \sum_{k=1}^{n} (Y_k-\mu)^2)</math>
 +
 +
 +
<math> = -\frac{n}{2} \frac{1}{\sigma^2}+\frac{1}{2} \sum_{k=1}^{n} (Y_k - \mu)^2(\hat{\sigma^2})^{-2} = 0 --- 2</math>
 +
 +
From Equation 1 and 2, we can find the ML estimator for the two parameters:
 +
 +
<center><math>\hat{\mu} = \frac{1}{n}\sum_{k=1}^n Y_k</math></center>
 +
<center><math>\hat{\sigma^2} = \frac{1}{n}\sum_{k=1}^n (Y_k - \hat{\mu})^2</math></center>
 +
 +
 +
 +
=Example 3: pD Gaussuan, <math>(\mu,R)</math> unknown=
 +
 +
Let <math>Y_1,...,Y_n</math> be i.i.d. Gaussian random vectors with distribution <math>N(\mu,R)</math> where <math>\mu \in \mathbb{R}^{p}</math> and <math>R \in \mathbb{R}^{p\times p}</math> is a symmetric positive definite matrix. Let <math>\theta = [\mu,R]</math> denote the parameter vector for the distribution.
 +
 +
 +
We would to know the density function for the entire observation first. In order to simplify the notation, we will represent the data as a matrix <math>Y</math>.
 +
 +
<center><math>Y = [Y_1,...,Y_n]</math></center>
 +
 +
Since the <math>n</math> observations, <math>Y_1,...,Y_n</math> are i.i.d., we know we can represent the density of <math>Y</math> as the following product of densities for each vector <math>Y_k</math>.
 +
 +
<center><math>p(x|\theta) = \prod_{k=1}^n \frac{1}{(2\pi)^{p/2}} |R|^{-1/2} \exp\{-\frac{1}{2} (y_k-\mu)^t R^{-1}(y_k-\mu)\}</math></center>
 +
 +
<center><math>p(x|\theta) = \frac{1}{(2\pi)^{pn/2}} |R|^{-\frac{n}{2}} \exp\{-\frac{1}{2}\sum_{k=1}^n (y_k-\mu)^t R^{-1}(y_k-\mu)\}</math></center>
 +
 +
 +
This expression may be furthur simplified by defining the following two new sample statistics,
 +
 +
<center><math> b = \sum_{k=1}^n y_k </math></center>
 +
<center><math> S = \sum_{k=1}^n y_ky_k^t = yy^t </math></center>
 +
 +
Where <math>yy^t</math> is the product of the matrix <math>y</math> with its transpose. Also,
 +
 +
<center><math> \sum_{k=1}^n (y_k-\mu)^t R^{-1}(y_k-\mu) </math></center>
 +
<center><math> = \sum_{k=1}^n (tr((y_k-\mu)^t R^{-1}(y_k-\mu))) </math></center>
 +
<center><math> = \sum_{k=1}^n (tr((y_k-\mu)^t (y_k-\mu)R^{-1})) </math></center>
 +
<center><math> = \sum_{k=1}^n tr((y_k y_k^t-2\mu y_k^t + \mu\mu^t)R^{-1}) </math></center>
 +
<center><math> = tr(\sum_{k=1}^ny_k y_k^t R^{-1}) - 2tr(\sum_{k=1}^ny_k^tR^{-1}\mu) + ntr(\mu^tR^{-1}\mu)</math></center>
 +
 +
where <math>tr\{A\}</math> denotes the trace of the matrix <math>A</math>, the summation of the diagonal elements.
 +
 +
In general, a statistic is any function of the observed data, but typically, statistics summarize the data in some useful way. Using these two statistics, the pdf of <math>Y</math> may be written as:
 +
 +
<center><math>p(x|\theta) = \frac{1}{(2\pi)^{pn/2}} |R|^{-\frac{n}{2}} \exp \{tr\{SR^{-1}\} +b^tR^{-1}\mu - \frac{n}{2}\mu^tR^{-1}\mu\}---3</math></center>
 +
 +
Maximizing this expression with respect to <math>\theta</math> then results in the ML estimate of the mean and covariance given by
 +
 +
<center><math> \arg\max_{[\mu, R]} p(y|\mu,R) = [\hat{\mu},\hat{R}]</math></center>
 +
 +
 +
The logarithm of the density function is given by:
 +
 +
<center><math> \log p(y|\theta) = -\frac{pn}{2}\log{2\pi} -\frac{n}{2}\log |R| - \frac{1}{2} tr(SR^{-1}) + b^t R^{-1}\mu - \frac{n}{2}\mu^tR^{-1}\mu</math></center>
 +
 +
Differentiating the log likelihood over <math>\mu</math> results in the following expressions.
 +
 +
<center><math> \bigtriangledown_{\mu} \log p_{\theta}(Y|\theta) =  b^t R^{-1} - \frac{n}{2}R^{-1}\hat{\mu}*2</math></center>
 +
 +
<center><math> \hat{\mu} = \frac{b}{n} = \frac{1}{n} \sum_{k=1}^n y_k</math></center>
 +
 +
However, finding the gradient of covariance matrix <math>R</math> is more complicated. Let
 +
 +
<center><math> M = \sum_{k=1}^n (y_k - \hat{\mu})(y_k-\hat{\mu})^t</math></center>
 +
<center><math> B = M^{\frac{1}{2}}R^{-1} M^{\frac{1}{2}}</math></center>
 +
 +
Then
 +
<center><math> R^{-1} = M^{\frac{1}{2}} B M^{-\frac{1}{2}}</math></center>
 +
<center><math> |R|^{-1} = |M^{-1}| |B|</math></center>
 +
 +
Using these two matrix, the density function can be written by:
 +
 +
<center><math> |R|^{-\frac{n}{2}} \exp \{-\frac{1}{2}tr(\sum_{k=1}^n(y_k-\hat{\mu})(y_k-\hat{\mu})^tR^{-1})\}</math></center>
 +
<center><math> = |M|^{-1/2}|B|^{n/2} \exp(-\frac{1}{2}tr(B))</math></center>
 +
 +
<center><math> \propto \prod_{i=1}^{p} \lambda_i ^{\frac{n}{2}}  \exp(-\frac{1}{2} \sum_{i=1}^{p} \lambda_i) = \prod_{i=1}^{p} \lambda_i ^{\frac{n}{2}}  \exp(-\frac{\lambda_i}{2}  )</math></center>
 +
 +
Therefore
 +
<center><math>\hat{R} = M^{1/2}B^{-1}M^{1/2} = \frac{1}{n} M</math></center>
 +
<center><math>\hat{R} = \frac{1}{n} \sum_{k=1}^{n} (y_k-\hat{\mu})(y_k-\hat{\mu})^t </math></center>
 +
<center><math>\hat{R} = \frac{1}{n} \sum_{k=1}^{n} y_ky_k^t + \sum_{k=1}^{n}\hat{\mu}{\hat{\mu}}^t - 2\sum_{k=1}^{n} y_k {\hat{\mu}}^t  </math></center>
 +
<center><math>\hat{R} = \frac{S}{n} + \hat{\mu}{\hat{\mu}}^t -\frac{2}{n}b{\hat{\mu}}^t = \frac{S}{n}- \hat{\mu}{\hat{\mu}}^t  </math></center>
 +
<center><math>\hat{R} = \frac{S}{n}- \hat{\mu}{\hat{\mu}}^t</math></center>
 +
 +
So from this we see that the ML estimate of the mean and covariance of a multivariate Gaussian distribution is simply given by the sample mean and covariance of the observations.
 
----
 
----
 
----
 
----

Latest revision as of 20:48, 30 April 2014


Maximum Likelihood Estimators and Examples
A slecture by Lu Zhang

Partially based on the ECE662 Spring 2014 lecture material of Prof. Mireille Boutin.




Outline of the slecture

  • Background and Introduction of ML estimator
  • Finding Maximum Likelihood Estimators (MLE)
  • Summary
  • References


Background and Introduction of ML estimator

Once we have decided on a model(Probability Distribution), our next step is often to estimate some information from the observed data. There are generally two parametric frameworks for estimating unknown information from data. We will refer to these two general frameworks as the Frequentist and Baysian approaches. One very widely used Frequentist estimator is known as the Maximum Likelihood estimator.


In the frequentist approach, one treats the unknown quantity as a deterministic, but unknown parameter vector, $ \theta \in \Omega $. So for example, after we observe the random vector $ Y \in \mathbb{R}^{n} $, then our objective is to use $ Y $ to estimate the unknown scalar or vector $ \theta $. In order to formulate this problem, we will assume that the vector $ Y $ has a probability density function given by $ p_{\theta}(y) $ where $ \theta $ parameterizes a family of density functions for $ Y $. We may then use this family of distributions to determine a function, $ T : \mathbb{R}^{n} \rightarrow \Omega $, that can be used to compute an estimate of the unknown parameter as

$ \hat{\theta} = T(Y) $

Notice, that since $ T(Y) $ is a function of random vector $ Y $, the estimate, $ \hat{\theta} $, is a random vector. The mean of the estimator, $ \bar{\theta} $, can be computed as

$ \bar{\theta} = E_{\theta}[\hat{\theta}] = \int_{\mathbb{R}^{n}} T(y)p_{\theta}(y)dy $

The difference between the mean of the estimator and the value of the parameter is known as the bias and is given by

$ bias_{\theta} = \bar{\theta} -\theta $


Similarly, the variance of the estimator is given by

$ var_{\theta} = E_{\theta}[(\hat{\theta} -\bar{\theta})^2] $

and it is easily shown that the mean squared error (MSE) of the estimate is then given by

$ MSE_{\theta} = E_{\theta}[(\hat{\theta}-\theta)^2] = var_{\theta} + (bias_{\theta})^2 $


Since the bias, variance, and the MSE of the estimator will depend on the specific value of $ \theta $, it is often unclear precisely how to compare the accuracy of different estimators. Even estimators that seem quite poor may produce small or zero error for certain values of $ \theta $. For example, consider the estimator which is fixed to the value $ \hat{\theta}=1 $, independent of the data. This would seem to be a very poor estimator, but it has an MSE of 0 when $ \theta=1 $.


An estimator is said to be consistent if for all $ \theta \in \Omega $, the MSE of the estimator goes to zero as the number of independent data samples, n, goes to infinity. If an estimator is not consistent, this means that even with arbitrarily large quantities of data, the estimate will not approach the true value of the parameter. Consistency would seem to be the least we would expect of an estimator, but we will later see that even some very intuitive estimators are not always consistent.


Ideally, it would be best if one could select an estimator which has uniformly low bias and variance for all values of $ \theta $. This is not always possible, but when it is we have names for such estimators. For example, $ \hat{\theta} $ is said to be an unbiased estimator if for all values of $ \theta $ the bias is zero, i.e. $ \theta = \bar{\theta} $. If in addition, for all values of $ \theta $, the variance of an estimator is less than or equal to that of all other unbiased estimators, then we say that the estimator is uniformly minimum variance unbiased(UMVU) estimator.


There are many excellent estimators that have been proposed through the years for many different types of problems. However, one very widely used Frequentist estimator is known as the maximum likelihood(ML) estimator given by

$ \hat{\theta} = \arg\max_{\theta \in \Omega} p_{\theta}(Y) $
$ = \arg\max_{\theta \in \Omega} \log p_{\theta}(Y) $

where the notation "argmax" denotes the value of the argument that achieves the global maximum of the function. Notice that these formulas for the ML estimate actually use the random variable $ Y $ as an argument to the density function $ p_{\theta}(y) $. This implies that $ \hat{\theta} $ is function of $ Y $, which in turn means that $ \hat{\theta} $ is a random variable.


When the density function, $ p_{\theta}(y) $, is a continuously differentiable function of $ \theta $, then a necessary condition when computing the ML estimate is that the gradient of the likelihood is zero.

$ \bigtriangledown_{\theta} p_{\theta}(Y)|_{\theta=\hat{\theta}} = 0 $


While the ML estimate is generally not unbiased, it does have a number of desirable properties. It can be shown that under relatively mild technical conditions, the ML estimate is both consistent and asymptotically efficient. This means that the ML estimate attains the Cramer-Rao bound asymptotically as the number of independent data samples,n, goes to infinity. Since this Cramer-Rao bound is a lower bound on the variance of an unbiased estimator, this means that the ML estimate is about as good as any unbiased estimator can be. So if one has no idea what values of $ \theta $ may occur, and needs to guarantee good performance for all cases, then the ML estimator is usually a good choice.




Finding the ML estimators

Example 1: 1D Gaussuan, $ \theta $ unknown, $ \sigma^2 $ known

Let $ \{Y_k\}_{k=1}^n $ be i.i.d. random variables with Gaussian distribution $ N(\theta,1) $ and unknown mean parameter $ \theta $. For this case, the logarithm of the density function is given by:

$ \log p_{\theta}(Y) = -\frac{1}{2\sigma^2}\sum_{k=1}^{n}(Y_k-\theta)^2 - \frac{n}{2} \log (2\pi\sigma^2) $


Local max are where gradient is zero. Differentiating the log likelihood results in the following expression.

$ \bigtriangledown_{\theta} \log p_{\theta}(Y) = \frac{d\log p_{\theta}(Y)}{d\theta} |_{\theta = \hat{\theta}} = \frac{1}{\sigma^2}\sum_{k=1}^{n}(Y_k - \hat{\theta}) = 0 $

From this we obtain the ML estimate for $ \theta $. So

$ \hat{\theta}=\frac{1}{n}\sum_{k=1}^{n} y_k $

Notice that, by the following argument, this particular ML estimate is unbiased.

$ E_{\theta}[\hat{\theta}] = E_{\theta}[\frac{1}{n}\sum_{k=1}^{n}Y_k] = \frac{1}{n}\sum_{k=1}^{n} E_{\theta}[Y_k] = \frac{1}{n} \sum_{k=1}^{n} \theta = \theta $

Moreover, for this special case, it can be shown that the ML estimator is also the UMVU estimator.



Example 2: 1D Gaussuan, $ (\mu,\sigma^2) $ unknown

Let $ \{Y_k\}_{k=1}^n $ be i.i.d. random variables with Gaussian distribution $ N(\mu,\sigma^2) $ and unknown mean parameter $ \theta $. For this case, the density function is giving by:

$ p_{\theta}(Y) = (2\pi\sigma^2)^{-\frac{n}{2}} \exp(- \frac{1}{2\sigma^2}\sum_{k=1}^{n}(Y_k - \mu)) $

The logarithm of the density function is given by:

$ \log p_{\theta}(Y) = -\frac{1}{2\sigma^2}\sum_{k=1}^{n}(Y_k-\mu)^2 - \frac{n}{2} \log (2\pi\sigma^2) $

Similarly, differentiating the log likelihood results in the following expressions.

$ \bigtriangledown_{\mu} \log p_{\theta}(Y) = \bigtriangledown_{\mu} (-\frac{n}{2}\log 2\pi-\frac{n}{2}\log \sigma^2 - \frac{1}{2\sigma^2} \sum_{k=1}^{n} (Y_k-\mu)^2) $

$ = -\frac{1}{2\sigma^2}\cdot 2\cdot (-1) \sum_{k=1}^{n}(Y_k-\hat{\mu}) $

$ = \frac{1}{\sigma^2} \sum_{k=1}^{n}(Y_k-\hat{\mu}) $

$ = \frac{1}{\sigma^2}(\sum_{k=1}^{n}Y_k -\hat{\mu}n) = 0 --- 1 $


$ \bigtriangledown_{\sigma^2} \log p_{\theta}(Y) = \bigtriangledown_{\sigma^2} (-\frac{n}{2}\log 2\pi-\frac{n}{2}\log \sigma^2 - \frac{1}{2\sigma^2} \sum_{k=1}^{n} (Y_k-\mu)^2) $


$ = -\frac{n}{2} \frac{1}{\sigma^2}+\frac{1}{2} \sum_{k=1}^{n} (Y_k - \mu)^2(\hat{\sigma^2})^{-2} = 0 --- 2 $

From Equation 1 and 2, we can find the ML estimator for the two parameters:

$ \hat{\mu} = \frac{1}{n}\sum_{k=1}^n Y_k $
$ \hat{\sigma^2} = \frac{1}{n}\sum_{k=1}^n (Y_k - \hat{\mu})^2 $


Example 3: pD Gaussuan, $ (\mu,R) $ unknown

Let $ Y_1,...,Y_n $ be i.i.d. Gaussian random vectors with distribution $ N(\mu,R) $ where $ \mu \in \mathbb{R}^{p} $ and $ R \in \mathbb{R}^{p\times p} $ is a symmetric positive definite matrix. Let $ \theta = [\mu,R] $ denote the parameter vector for the distribution.


We would to know the density function for the entire observation first. In order to simplify the notation, we will represent the data as a matrix $ Y $.

$ Y = [Y_1,...,Y_n] $

Since the $ n $ observations, $ Y_1,...,Y_n $ are i.i.d., we know we can represent the density of $ Y $ as the following product of densities for each vector $ Y_k $.

$ p(x|\theta) = \prod_{k=1}^n \frac{1}{(2\pi)^{p/2}} |R|^{-1/2} \exp\{-\frac{1}{2} (y_k-\mu)^t R^{-1}(y_k-\mu)\} $
$ p(x|\theta) = \frac{1}{(2\pi)^{pn/2}} |R|^{-\frac{n}{2}} \exp\{-\frac{1}{2}\sum_{k=1}^n (y_k-\mu)^t R^{-1}(y_k-\mu)\} $


This expression may be furthur simplified by defining the following two new sample statistics,

$ b = \sum_{k=1}^n y_k $
$ S = \sum_{k=1}^n y_ky_k^t = yy^t $

Where $ yy^t $ is the product of the matrix $ y $ with its transpose. Also,

$ \sum_{k=1}^n (y_k-\mu)^t R^{-1}(y_k-\mu) $
$ = \sum_{k=1}^n (tr((y_k-\mu)^t R^{-1}(y_k-\mu))) $
$ = \sum_{k=1}^n (tr((y_k-\mu)^t (y_k-\mu)R^{-1})) $
$ = \sum_{k=1}^n tr((y_k y_k^t-2\mu y_k^t + \mu\mu^t)R^{-1}) $
$ = tr(\sum_{k=1}^ny_k y_k^t R^{-1}) - 2tr(\sum_{k=1}^ny_k^tR^{-1}\mu) + ntr(\mu^tR^{-1}\mu) $

where $ tr\{A\} $ denotes the trace of the matrix $ A $, the summation of the diagonal elements.

In general, a statistic is any function of the observed data, but typically, statistics summarize the data in some useful way. Using these two statistics, the pdf of $ Y $ may be written as:

$ p(x|\theta) = \frac{1}{(2\pi)^{pn/2}} |R|^{-\frac{n}{2}} \exp \{tr\{SR^{-1}\} +b^tR^{-1}\mu - \frac{n}{2}\mu^tR^{-1}\mu\}---3 $

Maximizing this expression with respect to $ \theta $ then results in the ML estimate of the mean and covariance given by

$ \arg\max_{[\mu, R]} p(y|\mu,R) = [\hat{\mu},\hat{R}] $


The logarithm of the density function is given by:

$ \log p(y|\theta) = -\frac{pn}{2}\log{2\pi} -\frac{n}{2}\log |R| - \frac{1}{2} tr(SR^{-1}) + b^t R^{-1}\mu - \frac{n}{2}\mu^tR^{-1}\mu $

Differentiating the log likelihood over $ \mu $ results in the following expressions.

$ \bigtriangledown_{\mu} \log p_{\theta}(Y|\theta) = b^t R^{-1} - \frac{n}{2}R^{-1}\hat{\mu}*2 $
$ \hat{\mu} = \frac{b}{n} = \frac{1}{n} \sum_{k=1}^n y_k $

However, finding the gradient of covariance matrix $ R $ is more complicated. Let

$ M = \sum_{k=1}^n (y_k - \hat{\mu})(y_k-\hat{\mu})^t $
$ B = M^{\frac{1}{2}}R^{-1} M^{\frac{1}{2}} $

Then

$ R^{-1} = M^{\frac{1}{2}} B M^{-\frac{1}{2}} $
$ |R|^{-1} = |M^{-1}| |B| $

Using these two matrix, the density function can be written by:

$ |R|^{-\frac{n}{2}} \exp \{-\frac{1}{2}tr(\sum_{k=1}^n(y_k-\hat{\mu})(y_k-\hat{\mu})^tR^{-1})\} $
$ = |M|^{-1/2}|B|^{n/2} \exp(-\frac{1}{2}tr(B)) $
$ \propto \prod_{i=1}^{p} \lambda_i ^{\frac{n}{2}} \exp(-\frac{1}{2} \sum_{i=1}^{p} \lambda_i) = \prod_{i=1}^{p} \lambda_i ^{\frac{n}{2}} \exp(-\frac{\lambda_i}{2} ) $

Therefore

$ \hat{R} = M^{1/2}B^{-1}M^{1/2} = \frac{1}{n} M $
$ \hat{R} = \frac{1}{n} \sum_{k=1}^{n} (y_k-\hat{\mu})(y_k-\hat{\mu})^t $
$ \hat{R} = \frac{1}{n} \sum_{k=1}^{n} y_ky_k^t + \sum_{k=1}^{n}\hat{\mu}{\hat{\mu}}^t - 2\sum_{k=1}^{n} y_k {\hat{\mu}}^t $
$ \hat{R} = \frac{S}{n} + \hat{\mu}{\hat{\mu}}^t -\frac{2}{n}b{\hat{\mu}}^t = \frac{S}{n}- \hat{\mu}{\hat{\mu}}^t $
$ \hat{R} = \frac{S}{n}- \hat{\mu}{\hat{\mu}}^t $

So from this we see that the ML estimate of the mean and covariance of a multivariate Gaussian distribution is simply given by the sample mean and covariance of the observations.



Alumni Liaison

Prof. Math. Ohio State and Associate Dean
Outstanding Alumnus Purdue Math 2008

Jeff McNeal