\begin{center}
\textbf{ {2. MAXIMUM LIKELIHOOD ESTIMATION}}
\end{center}

{\textbf  2.1 Likelihood}

{ Likelihood} is one of the most fundamental concept in all types of
statistical inference.

{\textbf Definition 1} Suppose that $\bX$ has density function or
probability function $f(\bx ; \btheta)$. We have observed $\bX = \bx$.
Then the likelihood function with observation $\bx$ is defined as
\[
L(\btheta) \equiv
L(\btheta ; \bx) = f(\bx ; \btheta).
\]


{ \textbf{Density/probability function}: a function of $\bx$,
specifying the distribution of random variable $\bX$

\textbf{Likelihood}: a function of $\btheta$, reflecting information on
$\btheta$ contained in observation $\bx$ }


{\textbf Note}. A likelihood function represents the uncertainty on a
unknown non-random constant $\btheta$, and it is {\textbf not a density
or probability function}! It provides
\begin{quote}
  {\textbf  a rational degree of belief}, or\\
 {\textbf  an order of preferences}
\end{quote}
 on possible values of
the parameter $\btheta$. This can be seen more clearly in the simple
example on next slide.

In fact, a likelihood function is often defined up to a positive {
constant} --- the constant here refers to a quantity independent of
$\btheta$. But it may depending on $\bx$. (Note $\bx$ is a given
constant.)


{\textbf Example 1}. \label{bin} Suppose that $x$ is the number of
successes from a known number  $n$ of independent trials with unknown
probability of success $\pi$.  The probability function, and so the
likelihood function is
\[
L(\pi) =
f(x; \pi) = \begin{pmatrix}n\\x\end{pmatrix}\pi^x(1-\pi)^{n-x} .
\]
The likelihood function $L(\pi ; x)$ can be graphed as a function of
$\pi$. It changes shape for different values of $x$.  A likelihood
function for  $x=3$  when $n=10$ is shown in the Figure below.
\begin{center}\includegraphics[angle=-90,clip=T]{figure/fig21.eps}
\end{center}
 Notice that the likelihood function shown above  is {\em not} a
density function. It does not have an area of 1 below it.

We use the likelihood function to compare the plausibility of different
possible parameter values.  For instance,  the likelihood is much
larger for $\pi = 0.3$ than for $\pi = 0.8$, that is the data $x=3$
have a greater probability of being observed if $\pi = 0.3$ than if
$\pi = 0.8$.  This makes $\pi = 0.3$ much more {\textbf likely} as the
true value for $\pi$ than $0.8$.

{\textbf Note}. In the above argument, we do not need to calculate
exact probabilities under different values of $\btheta$. Only the order
of those quantities matters!


Let $X_1, \cdots, X_n$ be i.i.d. with pdf $f(\cdot, \btheta)$.
Write $\bX=(X_1, \cdots, X_n)^\tau$.
Then the likelihood function is
\[
L(\btheta) = L(\btheta; \bX) = \prod_{i=1}^n f(X_i, \btheta) ,
\]
which is a { product} of $n$ terms. Then the {\textbf log-likelihood
function} is
\[
l(\btheta) = l (\btheta; \bX) \equiv \log \{ L(\btheta; \bX)\}
= \sum_{i=1}^n \log \{ f(X_i, \btheta) \},
\]
which is a { sum} of $n$ terms.

This explains why log-likelihood functions are often used with
independent observations.


{\textbf Definition 2}. { \sl Likelihood Principle}

For two observed values $\bx$ and $\by$, if
\[
L(\btheta; \bx) \propto L(\btheta; \by),
\]
the inferences for $\btheta$ based on $\bx$ and $\by$ should
be the same.


{\textbf  2.2 Sufficiency and Data Compression}

If we start with $n$ observations and construct from them a
 $k$-dimensional statistic, if $k<n$ we expect that in general some
 information in the sample about $\btheta$ would be lost. The statistics
 is less complicated than the original sample because it is in a smaller
dimension. It reduces the data, and in general will reduce the
information about $\btheta$. A simple example may help to make clear what
is happening.



{\textbf Example 2} Based on a sample $\{ X_1 , \cdots , X_n\}$,
consider a sequence of statistics of increasing dimensions:
\begin{eqnarray*}
T_1(\bX)&=& \left(\sum_{i=1}^nX_i\right)\\
T_2(\bX)&=& \left(\sum_{i=1}^nX_i, \sum_{i=1}^n X_i^2\right)\\
T_3(\bX)&=& \left(\sum_{i=1}^nX_i, \sum_{i=1}^n X_i^2,\sum_{i=1}^nX_i^3\right)\\
&\dots&\\
T_{n-1}(\bX)&=& \left(\sum_{i=1}^nX_i, \sum_{i=1}^n
X_i^2,\sum_{i=1}^nX_i^3,\dots ,\sum_{i=1}^nX_i^{n-1}\right)\\
T_n(\bX)&=& \left(\sum_{i=1}^nX_i, \sum_{i=1}^n
X_i^2,\sum_{i=1}^nX_i^3,\dots ,\sum_{i=1}^nX_i^n\right)\\
\end{eqnarray*}
Here we see that $T_1$ is very simple and in 1 dimension, whereas $T_n$
is in $n$ dimensions.  Actually, from $T_n$ one can reconstruct the
original observations $X_1, X_2,\dots ,X_n$, except for their
ordering.  So  $T_n$ contains all the
information in the sample about $\btheta$.  But except in special cases
$T_1$ will not include all the information about $\btheta$ in the
sample.\\

If a statistic contains all the information in a sample about $\btheta$
we shall say that it is a \textbf{sufficient statistic} for $\btheta$.


{\textbf Example 3}
% \label{ex:suff}
If $X_1, X_2$ are independent observations from $N(\mu,1)$ then the
statistic $X_1+X_2$ contains all the information in the sample about
$\mu$.  This is intuitive, because knowing $X_1+X_2$ and $X_1-X_2$ is
equivalent to having both observations $X_1, X_2$, but $X_1-X_2$ has a
$N(0,2)$ distribution and so contains no information about $\mu$. This
conclusion is made more secure by the independence of $X_1+X_2$ and
$X_1-X_2$.\\



Now we make the idea of sufficiency more precise with one of the great
ideas from the work of R.A.Fisher.  We seek to define an idea of
`sufficiency' such that a statistic $T$ is `sufficient' for $\btheta$ if
it contains all the information about $\btheta$ that was in the whole
set of observations.


{\textbf Definition 3}.  {\sl Sufficient Statistic}\\
Suppose $\bX \sim f(\cdot, \btheta)$. $\bT(\bX)$ is said to be a
{\textbf sufficient statistic} for $\btheta$ if the conditional
distribution of the sample $\bX$ given $\bT(\bX)$ does not
depend on $\btheta$.\\

The likelihood principle implies the sufficiency principle which indicates
that we only need to use sufficient statistics in inference.

{\textbf Definition 4}.  { \sl Sufficiency Principle}\\
All sufficient statistics based on $\bX$ should lead to the same inferences
for $\btheta$.



{\textbf Theorem 1}. {\sl Factorisation Criterion}\\
Let $\bX \sim f(\bx, \btheta)$. Then $\bT=\bT(\bX) $ is a sufficient
statistic for $\btheta$ iff
\[
f(\bx;\btheta) = g(\bT(\bx),\btheta)h(\bx).
\]


\bigskip
\bigskip
\bigskip

We will not prove this result.  To see how a proof might go, note that
it is always possible, but not
necessary, to take $h(\bx)$ as the conditional density of $\bX$ given
$\bT(\bX)$ and $g(\bT(\bx),\btheta)$ as the density of
$\bT(\bX)$.


{\textbf Example 4}. Suppose that we have a random sample
$\bX=(X_1,\dots ,X_n)^\tau$ from $N(\mu,1)$.  Then the joint density is
\[
f_{\bX}(\bX;\btheta) = \frac{1}{(\sqrt{2\pi})^n}\exp
\left[-\frac{1}{2}\sum_{i=1}^n(X_i-\mu)^2\right]
\]
\[
= \exp\left[-\frac{1}{2}n(\bar{X}-\mu)^2\right] \frac{1}{(\sqrt{2\pi})^n}\exp
\left[-\frac{1}{2}\sum_{i=1}^n(X_i-\bar{X})^2\right].
\]
So we can take $T(\bX) = \bar{X}$ and $$g(\bar{X},\btheta)
 = \exp\left[-\frac{1}{2}n(\bar{X}-\mu)^2\right],$$ and $$h(\bX) =
 \frac{1}{(\sqrt{2\pi})^n}\exp\left[-\frac{1}{2}
 \sum_{i=1}^n(X_i-\bar{X})^2\right].$$

Two samples with the same sample mean should lead to exactly the same
inferences about $\mu$ if the sufficiency principle holds.


{\textbf Example 5}. Here is another example where a  sufficient
statistic is of the same dimension as the sample size, but is not the
whole sample. Let $\bX= (X_1, \cdots, X_n)$ be a sample from
$U(-\btheta,\btheta)$ and $\Theta = R^+$. Then the absolute values of
the observations $\mid X_1\mid, \mid X_2\mid,\dots ,\mid X_n\mid$ are
from $U(0,\btheta)$.   One loses the information on the signs of the
observations, but retains full information on $\btheta$.  Here the
statistic $\bT(X) =(\mid X_1\mid, \mid X_2\mid,\dots ,\mid X_n\mid)$ is
a sufficient statistic for $\btheta$.

In fact the conditional distribution
of $\bX$ given $\bT(\bX) = (T_1, T_2,\dots
,T_n)=(\mid X_1\mid, \mid X_2\mid,\dots ,\mid X_n\mid)$ is the
distribution where independently for all $i$, $X_i= T_i$ with
probability 0.5 and $X_i = -T_i$ with probability 0.5.  This
conditional distribution does not depend on $\btheta$, so the
definition is satisfied.  \\

{\textbf Example 6}. (Continuing Example) Another sufficient statistic
of $\btheta$ is $\max_{1\le i \le n} |X_i|$. To this end, note that
\[
f_{X_i}(x_i) = {\frac{1}{2} \btheta} I_{(-\btheta, \btheta) } (x_i)
=\left\{
\begin{array}{ll}
{1 \over 2 \btheta} \quad & |x_i|< \btheta,\\
0 & \mbox{otherwise}.
\end{array}
\right.
\]
Hence
\[
f_{\bX}(x_1, \cdots, x_n) = \prod_{i=1}^n f_{X_i}(x_i)
\]
\[
= \left\{
\begin{array}{ll}
{1 \over (2 \btheta)^n} \quad & \max_{1\le i\le n} |x_i| < \btheta,\\
0 & \mbox{otherwise}
\end{array}
\right.
\]
\[
= \; {1 \over (2 \btheta)^n} \, I ( \max_{1\le i\le n} |x_i| < \btheta).
\hspace{28mm}
\]
Now the conclusion follows from the factorisation criterion.

\bigskip

{\small How do you estimate $\btheta$?}


{\textbf Example 7}. Suppose that the $X_1, \cdots, X_n$ are iid $U(0,
\btheta)$, where $\btheta >0$. Then
\[
f(x_i; \btheta) = \left\{\begin{array}{ll}\frac{1}{\btheta} &x_i \in [0, \btheta]\\
                                          0&\text{ otherwise}
                         \end{array}
                 \right.
\]
which can be written
\[
 f(x_i; \btheta) =   \frac{1}{\btheta}\times I_{\{x_i\le\btheta\}}I_{\{x_i \ge 0\}}.
\]
The joint density is
\begin{align*}
f(\bx;\btheta)&=\frac{1}{\btheta^n}\prod_{i=1}^n   I_{\{x_i\le\btheta\}}  \prod_{i=1}^n I_{\{x_i \ge 0\}}\\
&=\frac{1}{\btheta^n}I_{\{\max_i x_i\le\btheta\}}I_{\{\min_i x_i \ge 0\}}\\
&=g\left(T(\bx);\btheta\right)\times h(\bx)
\end{align*}
where $T(\bx)= \max_i x_i$.  In this case the (at first sight) rather
surprising result is that $\max_i X_i$ is a  sufficient statistic for
$\btheta$.



{\textbf  2.3 Maximum likelihood estimator (MLE)}

The MLE is by far the most popular method for deriving estimators.

{\textbf Definition 5} --- {\textbf  MLE}\\
A {\em Maximum Likelihood Estimator} (MLE),
$\hat{\btheta}=\hat{\btheta}(\bX) \in \Theta$, of parameter $\btheta$ is an estimator
satisfying

\begin{quote}
 $L(\hat{\btheta};\bX)\ge L(\btheta ;\bX)$ for all $\btheta \in \Theta$, or
equivalently $l(\hat{\btheta};\bX)\ge l(\btheta ;\bX)$ for all $\btheta \in \Theta$.
\end{quote}

Obviously, a maximum likelihood estimator is  the most plausible value
for $\btheta$ as judged by the likelihood function.  In many cases
where $\Theta$ is continuous and the maximum does not occur at a
boundary of $\Theta$, $\hat{\btheta}$ {\textbf is often the solution}
of the equation
\[
s({\btheta}; \bX) = {\partial \over \partial \btheta} l(\btheta; \bX)
 =0.
\]
We call $s(\btheta) \equiv s({\btheta}; \bX)$ {\textbf  a score
function}.


{\textbf Example 8}. Suppose that $Y_1, Y_2,\ldots ,Y_n$ is a random
sample from $N(\mu, \sigma^2)$ where neither $\mu$ or $\sigma^2$ is
known.  Then we can find the maximum likelihood estimator from the
log-likelihood {\small
\begin{eqnarray*}
l(\mu ,\sigma^2 ) &=& -n\log\sqrt{2\pi}- n/2\log \sigma^2 -\sum_1^n(Y_i -\mu)^2/(2\sigma^2)\\
&
=& -n\log\sqrt{2\pi}- n/2\log \sigma^2-\sum_1^n(Y_i -\bar{Y})^2/(2\sigma^2)\\
&-&n(\bar{Y}-\mu)^2/(2\sigma^2).
\end{eqnarray*}}
This is maximised by choosing $\mu=\bar{Y}$, so $\hat{\mu}=\bar{Y}$ is
the MLE for $\mu$. It is easy to see
\[
E(\hat{\mu}) = {1\over n} \sum_{i=1}^n E(X_i) = \mu.
\]
Such a estimator is called {\textbf  unbiased}.

\askip

{\small {\textbf Note}. The above MLE for $\mu$ is equivalent to
minimise
\[
\sum_{i=1}^n(Y_i - \mu )^2.
\]
Therefore $\wh \mu$ is also the {\textbf  least squares estimator
(LSE)}. Note that LSE is derived from a simple empirical/geometric
rule, which makes \textbf{no use} of the underlying distribution. }




The {\textbf profile log-likelihood} remaining is
\[
l(\hat{\mu},\sigma^2 )=-n\log\sqrt{2\pi}+ (n/2)(\log \sigma^{-2}-
\wh \sigma^2 \sigma^{-2}),
\]
where $\wh \sigma^2 =\sum_1^n(Y_i -\bar{Y})^2/n$.
By the lemma below,
 the MLE for $\sigma^2$ is
$ \wh{\sigma}^2 .  $ Note that the MLE of $\sigma^2$ is {\em  biased}
since $$E(\hat{\sigma}^2) = (1-1/n) \sigma^2 \ne \sigma^2.$$

{\small {\textbf Lemma}. Define $L(x) = \log (x^{-1})  -b/x $, where
$n\ge 1$ and $b>0$ are constants. Then $L(b) \ge L(x)$ for all $x>0$. }

\bigskip
{\textbf Example 9}. Let $X_1, \cdots, X_n$ be i.i.d. Bernoulli($\pi$).
Then
\[
L(\pi) = \prod_{i=1}^n \pi^{X_i} (1-\pi)^{1-X_i} = \pi^{n\bar X} (1 -\pi) ^{n(1- \bar X)}.
\]
\[
l(\pi) = n\bar X \log \pi + n(1- \bar X) \log(1-\pi).
\]
Let $s(\pi) = {\partial \over \partial \pi} l(\pi) = 0$, leading to $\wh \pi = \bar X$.



{\textbf Example 10}. Suppose that  $Y_1, Y_2,\ldots ,Y_n$ is a random
sample from an exponential distribution with density function
$e^{-(y-\theta)}$ for $y\ge\theta$.  This is the usual exponential
distribution shifted to start at $\theta$.  The Likelihood is
\[
L(\theta ; \bY) = e^{-n(\bar{Y}-\theta)}I_{\{(\theta ,\infty)\}}(Y_{(1)}),
\]
where $Y_{(1)}$ is the smallest observation.  This likelihood is zero
for $\theta > Y_{(1)}$ and increases in $\theta$ for $\theta \le
Y_{(1)}$.  So the MLE $\hat{\theta}=Y_{(1)}$, which is a boundary
maximum.\\

\centerline{\includegraphics{figure/fig22.ps}}


{ Invariance property of MLEs}

Suppose $\bX \sim f(\bx, \btheta)$, and $\bmu= \bmu(\btheta)$ is a
1-1 transform. Let $\wh \btheta$ be the MLE for $\btheta$, i.e.
\[
\wh \btheta = \arg\max_{\btheta} f(\bX, \btheta).
\]
It is obvious
to see that the MLE for $\bmu$ is
\[
\wh \bmu = \bmu(\wh \btheta).
\]



{\textbf  2.4 Numerical computation of MLEs}

In modern statistical applications, it is  typically difficult to find explicit
analytic forms for the maximum likelihood estimators. These estimators are
found more often by iterative procedures built into computer software. An
iterative scheme starts with some guess at the MLE and then steadily
improves it with each iteration.  The estimator is considered found when it
has become numerically  stable. Sometimes the iterative procedures become
trapped at a local maximum which
is not a global maximum.  There may be a very large number of parameters in
a model, which makes such local entrapment much more common.


\bigskip

{\textbf  Newton-Raphson Scheme}

Suppose that the log-likelihood function function $l(\btheta)$ is
sufficiently smooth. Then
\[
s(\wh \btheta) =0,
\]
where $\wh \btheta$ is the MLE and $s(\btheta) = {\partial \over \partial \theta}
 l(\btheta)$ is the
score function. Let
\[
\dot s(\btheta) = \ddot  l(\btheta) = {\partial^2 \over \partial \theta
\partial \theta^\tau} l(\btheta).
\]

Suppose $\wh \btheta$ is close to the true value $\btheta_0$.
By a simple Taylor expansion,
\[
s(\btheta_0) =  \dot s(\btheta_0) (\btheta_0  - \wh \btheta) + o_p(||\wh
\btheta - \btheta_0||).
\]
This leads to the approximation
$$
\wh \btheta  \approx \btheta_0 - \{ \dot s( \btheta_0) \}^{-1} s(\btheta_0).
$$
Since $\theta_0$ is unknown, we use iterative estimators {
\begin{equation} \label{nr}
\wh \btheta_{k+1} = \wh \btheta_k - \{ \dot s( \wh \btheta_k) \}^{-1} s(\wh
\btheta_k)
\end{equation}
}
for $k=1, 2, \cdots$, where $\wh \btheta_0$ is a prescribed initial value.
We define $\wh \btheta = \wh \btheta_j$ if $\wh \btheta_j$ and
$\wh \btheta_{j-1}$ differ by a small amount.

\bigskip

{\textbf  Fisher Score method}: replace $\dot s( \wh \btheta_k)$ in
(\ref{nr}) by $ E_{\btheta} \{ \dot s( \btheta ) \}$ under $\btheta =
\wh \btheta_k$.



Like most iterative algorithms, the choice of appropriate {\textbf
 initial values is important} to ensure the convergence to
right limits. In practice multiple initial values are often used.

The differences between the Newton-Raphson and score methods are subtle.
We make observations below
\begin{itemize}
\item
The convergence of the Newton-Raphson algorithm is often faster
when both algorithms converge

\item
The radius of convergence for the score method is often larger, making
the choice of initial values less important for the score method.
\end{itemize}





{\textbf  {2.5 EM  algorithms}}

{\textbf  Goal}: to find the MLE $\wh \btheta = \wh \btheta(\bY)$ for
$\btheta$ from the likelihood based on data $\bY$:
\[
L(\btheta; \bY) = f_\bY (\bY, \btheta),
\]
while the `{\sl  complete}' data $\bX^\tau = (\bY^\tau, \bZ^\tau)$
contain a `{\sl  missing}' component $\bZ$. The likelihood based on the
 complete data is
\[
L(\btheta; \bX) = f_\bX (\bX, \btheta).
\]

\bigskip

{{\textbf EM} ({\sl E}xpectation and {\sl M}aximisation) {\textbf
algorithms}
 }

\begin{quote}
$\bullet$ {\sl  E-step}: compute the conditional expectation
\[
Q(\btheta) = Q(\btheta | \bY, \btheta_0) \equiv E \{ \log L(\btheta; \bX) | \bY,
\btheta_0\}
\]
$\bullet$ {\sl  M-step}: maximise $Q(\btheta)$ to give an updated value
$\btheta_1$
\end{quote}

then go to the E-step using $\btheta_0= \btheta_1$, and keep iterating until
convergence.
The limit of $\btheta_0$ is taken as $\wh \btheta(\bY)$.




{\textbf Example 11}. The genetic example from (Rao 1973, p.396)
assumes that the phenotype data
\begin{eqnarray*}
\bY &=& (Y_1, Y_2, Y_3, Y_4)^\tau \; \sim \\
 && \mbox{Multinomial}\big(4; {1\over 2}
+ {\theta \over 4},
{1-\theta \over 4}, {1-\theta \over 4}, {\theta \over 4} \big),
\end{eqnarray*}
where $\theta \in (0,1)$.
Then log-likelihood is
{\small
\[
l(\theta, \bY) = Y_1 \log(2+\theta) + (Y_2+Y_3) \log(1-\theta) + Y_4 \log\theta + C,
\]}
which {\textbf  does not yield a closed form $\wh \theta$}.

Now we treat $\bY$ as incomplete data from $\bX=(X_1, \dots, X_5)^\tau$ with
multinomial probabilities
\[
\big( {1\over 2},\; {\theta\over 4},\; {1-\theta \over 4},\; {1-\theta
\over 4}, \;{\theta \over 4} \big).
\]
Then
\[
Y_1 = X_1 + X_2, \quad Y_i = X_{i+1} \;\;\; {\rm for} \; i=2,3,4.
\]
The log-likelihood of based on $\bX$ is
{\small
\[
l(\theta, \bX) = (X_2+X_5) \log \theta + (X_3 + X_4)\log(1-\theta) + C,
\]}
which {\textbf  readily yields}
\[
\wh \theta (\bX) = { X_2 + X_5 \over X_2 + X_3 + X_4 + X_5}.
\]

Now the E-step is to find
{\small
\begin{eqnarray*}
Q(\theta) &=&  E\{ l(\theta, \bX) | \bY, \theta_0\}\\
&=& \log \theta E( X_2 + X_5 | \bY, \theta_0) +
\log (1- \theta) E( X_3 + X_4| \bY, \theta_0)\\
&=& (\wh X_2 + Y_4) \log \theta
+ (Y_2 + Y_3 ) \log (1- \theta),
\end{eqnarray*}
}
where
$\wh X_2 = E(X_2 | \bY, \theta_0).$
Since the conditional distribution of $X_2$ given $Y_1 (= X_1 + X_2)$ is
a binomial distribution with $n=Y_1$ and
\[
p= {\theta_0/4 \over 1/2 + \theta_0/4} = { \theta_0 \over 2 + \theta_0}.
\]
Hence
\begin{equation} \label{ex111}
\wh X_2 = n p = {Y _1 \theta_0 \over 2 + \theta_0}.
\end{equation}

The M-step leads to
\begin{equation} \label{ex112}
\theta_1 = { \wh X_2 + Y_4 \over \wh X_2 + Y_4  + Y_2 + Y_3 }.
\end{equation}

\bigskip

{\small For $Y=(125, 18, 20, 34)$ and the \textbf{initial value}
$$ \theta_0 = 4\times 34 /(125+18+20+34)$$
which is a relative frequency estimate, the first 5 iterations
between (\ref{ex111}) and (\ref{ex112}) are
0.690, 0.635, 0.628, 0.627 and 0.627, giving the MLE
$\wh \theta = 0.627$.
}


What can be said about {\sl convergence properties} of the EM algorithm?

Let $\btheta_0$ be an arbitrary initial value and $\btheta_1$ be the updated
value obtained from applying the iteration once.
Then it can be shown that
\[
L( \btheta_1; \bY) \ge L( \btheta_0; \bY).
\]
Unfortunately, it does not imply that the iterations will always lead to
the MLE eventually.

It is important to choose appropriate initial values to ensure the
algorithm converges to the MLE. (This is also true for both
Gaussian-Raphson and score methods!) In practice, it is a good idea to
use {\textbf  a variety of initial values}.


Further discussion on the convergence of EM algorithms, see
Wu (1983) {\sl Annals of Statistics, Vol.11, pp.95-103} and
\S 12.4 of Pawitan~(2001).




{\textbf  General comments on EM algorithms}:

{\small
\begin{itemize}
\item
The EM algorithm is a general procedure for
computing MLEs. It is not really a numerical algorithm. The calculation
of M-Step typically involves other numerical algorithms such as
Newton-Raphson and the score methods.

\item
It can be applied when some data are {\em genuinely} missing. It can
also be applied when the missing information is  merely a concept based
on which we transform a difficult optimisation problem into a sequence
of easier problems; see Example~11 above. This is particularly relevant
when $\wh \btheta(\bY)$ is difficult to calculate while $\wh \btheta(\bX)$
is easier to obtain.

\item
The convergence of EM algorithm may be very slow, depending on the amount
of missing information. It is often
slower than the Newton-Raphson algorithm.
\end{itemize}
}





{\textbf Example 12}. { Mixture distributions}\\
Let $Y_1, \cdots, Y_n$ be i.i.d. with a mixture pdf
\[
f(y) = \sum_{j=1}^k \alpha_j f_j(y, \blambda_j),
\]
where $f_j$ are pdfs, $\alpha_j\ge 0$ and $\sum_j \alpha_j =1$.
The parameter $\btheta$ contains all $\alpha_j$ and $\blambda_j$.
This model becomes important because
\begin{quote}
1. it represents heterogeneous data well since each $f_j$ represents
one heterogeneous component, and\\
2. it provides very good approximations to a large class of distributions.
\end{quote}


The likelihood based on $\bY=(Y_1, \cdots, Y_n)^\tau$ is
\[
L(\btheta; \bY) = \prod_{i=1}^n\Big\{ \sum_{j=1}^k \alpha_j f_j(Y_j, \blambda_j)
\Big\}.
\]
Maximising this likelihood is difficult, due to the presence of the
summations, which reflect the fact that we are typically lacking in
knowledge of which component any particular sample value comes from.
This is the missing information!

Let $\bX_i^\tau = (Y_i, \bZ_i^\tau)$, where $\bZ_i= (Z_{i1}, \cdots, Z_{ik})^\tau$
is a $k\times 1$ vector with
a 1 in the position corresponding to the component of the mixture that $Y_i$ comes
from, and 0 elsewhere.

Let
$\be_j $ be the $k\times 1$ vector with the $j$-th component 1 and all
the other components 0.
The joint probability-density function for $(\bZ_1, Y_1)$ is
\begin{eqnarray*}
&& P\{ Z_1 = \be_\ell,\; Y_1 \in [y_1, y_1 + dy_1) \}/dy_1\\
&=& P (Z_1 = \be_\ell) \, P \{ Y_1 \in [y_1, y_1 + dy_1)|Z_1 = \be_\ell\} /dy\\
&=& \alpha_\ell f_\ell (y_1, \blambda_\ell)
\; = \; \prod_{j=1}^k \alpha_j^{e_{\ell j}} f_j(y_1, \blambda_j) ^{e_{\ell j}},
\end{eqnarray*}
where $\be_\ell = (e_{\ell 1}, \cdots, e_{\ell k} )^\tau$.

 Let $\bX^\tau = (\bX_1^\tau, \cdots, \bX_n^\tau)$.
\[
L(\btheta; \bX) = \prod_{i=1}^n \prod_{j=1}^k \alpha^{Z_{ij}}_j f_j (Y_i,
\blambda_j)^{Z_{ij}},
\]
{\small
\[
l(\btheta; \bX) = \sum_{i=1}^n
\bZ_i^\tau
\left\{
\left(
\begin{array}{c}
\log \alpha_1\\
\vdots\\
\log \alpha_k
\end{array}
\right)
+
\left(
\begin{array}{c}
\log f_1(Y_i, \blambda_1)\\
\vdots\\
\log f_k(Y_i, \blambda_k)
\end{array}
\right)
\right\}.
\]
}
Note that
{\small
\[
E(\bZ_i | \bY, \btheta_0 )
=
\Big( {
\alpha_1^0 f_1(Y_i | \blambda_1^0) \over \sum_{\ell = 1}^k
\alpha_\ell^0 f_\ell(Y_i | \blambda_\ell^0)} ,
\cdots,
{\alpha_k^0  f_k(Y_i | \blambda_k^0)\over \sum_{\ell = 1}^k
\alpha_\ell^0 f_\ell(Y_i | \blambda_\ell^0)}
\Big)^\tau ,
\]}
which are constants as far as the M-step is concerned.


Now the E-step implies that
{\small
\[
Q(\btheta) = \sum_{i=1}^n
E(\bZ_i^\tau | \bY, \btheta_0 )
\left\{
\left(
\begin{array}{c}
\log \alpha_1\\
\vdots\\
\log \alpha_k
\end{array}
\right)
+
\left(
\begin{array}{c}
\log f_1(Y_i, \blambda_1)\\
\vdots\\
\log f_k(Y_i, \blambda_k)
\end{array}
\right)
\right\},
\]
}

and the M-step requires to maximise this function, which is much
easier than minimizing $l(\theta, \bY)$ directly for, for example,
normal pdf $f_j$'s.




{\textbf  2.6 Evaluating estimation}

To measure the accuracy of an MLE or, more general,  any estimation
procedure, we need to define some measures for the goodness (or badness)
of an estimator.

Let $\wh \btheta = \wh \btheta(\bX)$ be an estimator of $\btheta$, and
$\btheta_o$ be the (unknown) {\sl true value} of $\btheta$. Note that
 {
\begin{quote}
(i) exact estimation error $\wh \btheta  - \btheta_o$ is unknown, and\\
(ii) $\wh \btheta$ is a random variable
\end{quote}}
 we have to gauge the error  {
\begin{quote}
(i)
in terms of a probability average, and\\
(ii)
for all possible values of $\btheta_o \in \Theta$.
\end{quote}}

Let $P_\btheta$, $E_\btheta$ and $\var_\btheta$ denote the probability distribution,
expectation and variance under $\btheta_o = \btheta$.


{\textbf  Bias}: Bias$_\btheta( \wh \btheta)= E_\btheta(\wh \btheta) -
\btheta$

{\textbf  Variance}: $ \var_\btheta(\wh \btheta) $

{\textbf  Standard error}: $\{ \var_\btheta(\wh \btheta) \}^{1/2}$

{\textbf  Mean square error (MSE)}: $ E_\btheta ( \wh \btheta -
\btheta)^2$

{\textbf  Mean absolute error (MAE)}: $ E_\btheta|\wh \btheta -
\btheta|$

\bigskip

Note that
\begin{itemize}
\item
standard error is a meaningful measure of accuracy
 for unbiased estimators only, and
\item MSE (it its squared-root) should be used in general as {
\[
\MSE_\btheta (\wh \btheta) = \{\bias_\btheta(\wh \btheta)\}^2 + \var_\btheta(\wh \btheta).
\]}
\end{itemize}

Ideally we would seek for the estimator which minimises MSE or MAE
{\textbf for all $\btheta \in \Theta$}  over all possible candidate
estimators. Unfortunately such a global optimum rarely exists. However
if we confine to some subclass of estimators, the MLE is often optimal
or asymptotically optimal.

The MSE is most frequently used largely due to is { technical
tractability} while the MAE leads to estimators which is { more robust}
against outliers in observations.






{\textbf  Fisher Information}

Suppose $\bX \sim f(\bx, \btheta)$. The score function is
\[
s(\btheta) = \dot l (\btheta; \bX) = {\partial \over \partial \btheta}
\log \{ f(\bX, \btheta)\}.
\]
We assume certain regularity conditions so that we can take derivatives
under the integral sign.

\textbf{Mean of $s(\btheta)$}:
\begin{eqnarray*}
E_\btheta\{ s(\btheta)\} &=& \int s(\btheta) f(\bx, \btheta) d \bx\\
&=& \int {\partial \over \partial \btheta}
\log \{ f(\bx, \btheta)\} f(\bx, \btheta) d \bx\\
&=& \int {\partial \over \partial \btheta} f(\bx, \btheta) d \bx\\
&=& {\partial \over \partial \btheta} \int f(\bx, \btheta) d \bx \; = \; 0.
\end{eqnarray*}



\textbf{Variance of $s(\btheta)$} --- {\textbf  Fisher information}
matrix:
\begin{eqnarray*}
\calI(\btheta) & = & \calI_\bX(\btheta) = \var_\btheta\{ s(\btheta)\} \\[1ex]
&=& E _\btheta\{ s(\btheta) s(\btheta)^\tau \} \\
&=&
- E _\btheta\big[ {\partial^2 \over \partial \btheta \partial \btheta^\tau }
\log \{ f(\bX, \btheta)\} \big],
\end{eqnarray*}
because
\begin{eqnarray*}
&&
E _\btheta\big\{ {\partial^2 \over \partial \btheta \partial \btheta^\tau }
l(\btheta) \big\} \\[1ex]
&=&
\int
{ \ddot L(\btheta) L (\btheta) - \dot L(\btheta) \dot L(\btheta)^\tau
\over  L (\theta) } d \bx \\[1ex]
&=& \int \ddot L(\btheta) d \bx - \int { \dot L(\btheta) \dot L(\btheta)^\tau
\over L (\theta) } d \bx\\[1ex]
&=& - \int { \dot L(\btheta) \dot L(\btheta)^\tau
\over L (\theta) } d \bx\\[1ex]
&=& - \int s(\btheta) s(\btheta)^\tau f(\bx, \btheta) d \bx \\[1ex]
&=&  - E _\btheta\{ s(\btheta) s(\btheta)^\tau \}.
\end{eqnarray*}

Fisher information $\calI(\btheta)$ measures { \textbf the information
on $\btheta$ contained in data $\bX$}. Further if $\bX = (X_1, \cdots,
X_n)^\tau$, and $X_1, \cdots, X_n$ are independent,
\[
\calI(\btheta) = \calI_\bX(\btheta) = \sum_{j=1}^n \calI_{X_j} (\btheta).
\]


\askip

For $\btheta = \theta$ is a scalar, the Fisher information is
\[
\calI(\theta) = E_\theta \{ s(\theta)^2 \} = - E_\theta \{ \ddot l (\theta) \}.
\]


\askip

{\textbf Theorem 2}. {\em (Cram\'er-Rao inequality)}\\
Let $\bX \sim f(\cdot, \theta)$ which satisfying some regularity conditions.
Let $T=T(\bX)$ be a statistic with $g(\theta) = E_\theta (T)$. Then
for any $\theta \in \Theta$,
\[
\var_\theta (T) \ge \{ \dot g(\theta) \}^2/\calI(\theta).
\]



The Cram\'er-Rao inequality specifies {\textbf a lower bound} for any
{\sl unbiased estimator} for the parameter $g(\theta)$. When the
equality holds, $T$ is the {\textbf  minimum variance unbiased
estimator (MVUE)} of $g(\theta)$.


An MVUE is a member in the class of all unbiased estimators with the
minimum variance. {\sl  Not all MVUEs can attain the Cram\'er-Rao lower
bounds}.




{\textbf Example 13}. Let $X_1, \cdots , X_n$ be a sample from $N(\mu ,
\sigma^2)$.  We consider estimators for $\mu$, treating $\sigma^2$ as
known. The score function (for one observation) is
\begin{eqnarray*}
&&
s(\mu ; X_1) = \frac{\partial}{\partial\mu } \log [e^{-\frac{1}{2\sigma^2}(X_1-\mu)^2}/\sqrt{2\pi\sigma^2}]\\
&=&\frac{\partial}{\partial\mu }[-\frac{1}{2\sigma^2}(X_1-\mu)^2]
=(X_1-\mu)/\sigma^2.
\end{eqnarray*}
Note $\ddot l(\mu) = \dot s(\mu) = -\sigma^{-2}$.
Hence the Fisher information based on a single observation is
$\calI_{X_1} (\mu) = \sigma^{-2}$. Therefore
\[
\calI(\mu) =
\calI_{X_1, \cdots, X_n}(\mu) = n/\sigma^2.
\]
For any unbiased estimator $\wh \mu$ for $\mu$, it holds that
\[
\var_\mu (\wh \mu) \ge \sigma^2/n,
\]
which is the variance of $\bar X$. Hence $\bar X $ is the MVUE for $\mu$.


{\textbf  Asymptotic properties of MLEs}

Let $X_1, \cdots , X_n$ be i.i.d. with pdf $f(\cdot, \btheta)$. Write
\[
l(\btheta) = l(\btheta; \bX) = \sum_{j=1}^n \log f(X_j, \btheta).
\]
Let $\wh \btheta $ be the MLE which maximises $l(\btheta)$.
Suppose $f$ fulfils certain regularity conditions.

(a) \textbf{ Consistency}.\\
The MLE is consistent in the sense that as $n\to \infty$,
\[
P_\btheta \{ || \wh \btheta - \btheta || > \ve \} \to 0
\]
for any $\ve>0$.

Consistency requires that an estimator converges to the parameter to be
estimated. It is a very mild and modest condition that any reasonable estimator
should fulfil. The consistency condition is often used to {\sl rule out
bad estimators}.



(b) \textbf{ Asymptotic normality}\\
As $n \to \infty$,
\[
n^{1/2} ( \wh \btheta - \btheta ) \wcon N\big(0, \; \{\calI_{X_1}(\btheta)\}^{-1} \big).
\]

For large $n$, it holds {\textbf approximately} that
\[
\wh \btheta \sim N \big( \btheta,\; \{\calI_{X_1}(\btheta)\}^{-1}/n \big).
\]
Therefore {\sl  asymptotically} the MLE is unbiased and attains the
Cram\'er-Rao lower bound. Any estimator fulfilling this condition is
called {\textbf  efficient}.

{\textbf  An approximate standard error} of the $j$-th component of
$\wh \btheta$ is the square-root of the $(j,j)$-th element of
$\{\calI_{X_1}(\btheta)\}^{-1}$ divided by $n^{1/2}$.





{\textbf  {2.7 Bootstrapping MSEs }}

An MSE provides a measure for the accuracy of the estimator.
The standard errors of an MLE derived from the asymptotic normality typically depend on
unknown parameters, which are replaced by their estimators in practice.
Alternatively we may estimate the MSE by a simulation method
called Bootstrap.

Let $X_1, \cdots , X_n$ be i.i.d. with pdf $f(\cdot, \theta)$. Let
$$\wh \theta = T(X_1, \cdots, X_n)$$ be an estimator. The
{\sl goal} here is to estimate
\[
\nu \equiv \{ \MSE_{\theta_o}(\wh \theta)\}^{1/2},
\]
where $\theta_o$ is the true value.

If we {\sl knew} $f(\cdot, \theta_o)$ completely, $\nu$ is known in
principle, and  may be estimated
easily via a repeated sampling as follows.   We draw $B$ independent samples of
size $n$ from  $f(\cdot, \theta_o)$. For each sample, we calculate $\wh
\theta$, obtaining $\wh \theta_1, \cdots, \wh \theta_B$.
 Then the sample root-MSE
\[
\big\{ {1 \over B} \sum_{b=1}^B ( \wh \theta_b - \theta_o)^2 \big\}^{1/2}
\]
is a reasonable estimator for
$\nu$. By the LLN,  this estimator converges to $\nu$
as $B\to \infty$.

The {\textbf basic idea of bootstrap} is to adopt the above sampling
procedure in the so-called {\sl  bootstrap world}: now the population
is $f(\cdot, \wh \theta)$ which is \textbf{known}. We draw a sample
denoted as $(X_1^*, \cdots,\\ X_n^*)$ from this distribution. Define
the {\sl bootstrap version} of the estimator
\[
\wh \theta^* = T(X_1^*, \cdots, X_n^*).
\]
Then the quantity
\[
\nu^* =  \{\MSE_{\wh \theta} ( \wh \theta^* )\}^{1/2}
\]
is known in principle since the distribution $f(\cdot, \wh \theta)$ is
completely known. Define $\nu^*$ as a {\textbf bootstrap estimator} for
$\nu$.

In practice, we draw $B$ sets samples from $f(\cdot, \wh \theta)$, forming
$B$ bootstrap versions of estimator $\wh \theta^*_1, \cdots ,
\wh \theta^*_B$. The $\nu^*$ is calculated as
\[
\nu^* =\big\{ {1 \over B} \sum_{j=1}^B ( \wh \theta^*_j -
\wh \theta)^2 \big\}^{1/2}.
\]

\askip
\askip
\askip
\askip
\askip
\askip
\askip

\begin{center}
{

\psset{framearc=.4}

\begin{psmatrix}[rowsep=30pt,colsep=68pt]
Real World & { Bootstrap World}\\
\psframebox{
$
\begin{array}{c}
{\rm Population}\\
f(\cdot, \theta_o)\\
{\rm unknown}
\end{array}
$} &
\psframebox{
\
$
\begin{array}{c}
{\rm Population}\\
f(\cdot, \wh\theta)\\
{\rm known}
\end{array}
$ }\\
\psframebox{
$
\begin{array}{c}
{\rm Estimator}\\
\wh\theta= T(X_1, \cdots, X_n)
\end{array}
$} &
\psframebox{
\
$
\begin{array}{c}
{\rm Estimator}\\
\wh\theta^*= T(X_1^*, \cdots, X_n^*)
\end{array}
$ }\\
\psframebox{
$
\begin{array}{c}
{\rm Error}\\
\wh \theta - \theta_o\\
{\rm unknown}
\end{array}
$} &
\psframebox{

$
\begin{array}{c}
{\rm Error}\\
\wh \theta^* - \wh \theta\\
{\rm known}
\end{array}
$ } \psset{linestyle=solid} \ncline{2,1}{3,1} \ncline{3,1}{4,1}
\ncline{2,2}{3,2} \ncline{3,2}{4,2}
\ncline{<-}{4,1}{4,2}\naput{\tiny{Estimating}}
\end{psmatrix}
}
\end{center}


{ Bootstrap is a powerful tool for statistical inference}. It has
different forms for different applications. The one introduced here is
in the form of {\sl parametric bootstrap}. The special monographs on
this topic include:

{\small

\begin{itemize}
\item
Hall, P. (1992). The bootstrap and Edgeworth expansion. Springer, New York.

\item
Efron, B. and Tibshirani, R.J. (1993). An introduction to the bootstrap.
Chapman and Hall, New York.

\item
Shao, J. and Tu, D. (1995). The jackknife and bootstrap. Springer, New York.

\item
Davison, A.C. and Hinkley, D.V. (1997). Bootstrap methods and their application.
Cambridge University Press, Cambridge.
\end{itemize}
}
