\begin{center}
{\bf \colorbox{grey}{2. MAXIMUM LIKELIHOOD ESTIMATION}}
\end{center}

{\bf \color{dblue} 2.1 Likelihood}

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

{\bf 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).
\]


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

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


{\bf Note}. A likelihood function 
represents the uncertainty on a unknown non-random constant $\btheta$, and 
it is {\bf not a density or probability function}! 
It provides
\vspace{-4mm}
\begin{quote}
  {\bf \color{red} a rational degree of belief}, or\\
 {\bf \color{red} 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
{\color{dgreen} constant} --- the constant here refers to a quantity
independent of $\btheta$. But it may depending on $\bx$. (Note $\bx$ is
a given constant.)

\newpage
{\bf 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 a $x=3$  when $n=10$ is shown in the Figure below.
\centerline{\epsfig{figure=figure/fig21.ps, angle=-90}}
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 {\bf likely} as
the true value for $\pi$ than $0.8$.

{\bf 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!

\newpage
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 {\color{blue} product} of $n$ terms.
Then the {\bf 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 {\color{blue} sum} of $n$ terms.

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

{\bf Definition 2}. {\color{dblue} \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.

\newpage
{\bf \color{dblue} 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.

 
\newpage
{\bf 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 \emph{sufficient statistic} for $\btheta$.


{\bf 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$.\\


\newpage
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.


{\bf Definition 3}.  {\color{dblue}\sl Sufficient Statistic}\\
Suppose $\bX \sim f(\cdot, \btheta)$.
$\bT(\bX)$ is said to be a {\bf 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.

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

\newpage

{\bf Theorem 1}. {\color{dblue}\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)$.

\newpage
{\bf 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.

\newpage
{\bf 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.  \\

{\bf 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) = {1 \over 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$?}

\newpage
{\bf 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$.

\newpage

{\bf \color{dblue} 2.3 Maximum likelihood estimator (MLE)}

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

{\bf Definition 5} --- {\bf \color{dblue} MLE}\\
A {\em Maximum Likelihood Estimator} (MLE),
$\hat{\btheta}=\hat{\btheta}(\bX) \in \Theta$, of parameter $\btheta$ is an estimator
satisfying
\vspace{-4mm}
\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}$ {\bf 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)$ {\bf \color{red} a score function}.

\newpage
{\bf 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 {\bf \color{dgreen} unbiased}.

\askip

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


\newpage

The {\bf \color{blue}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 \color{dgreen} biased} 
since $$E(\hat{\sigma}^2) = (1-1/n) \sigma^2 \ne \sigma^2.$$

{\small
{\bf 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
{\bf 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$.

\newpage

{\bf 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{\epsfig{figure=figure/fig22.ps}}

\newpage
{\color{dmagenta} 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).
\]

\newpage

{\bf \color{dblue} 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

{\bf \color{dgreen} 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
{\color{blue}
\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

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

\newpage

Like most iterative algorithms, the choice of appropriate {\bf
\color{dred} 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}



\newpage

{\bf \fcolorbox{blue}{yellow} {2.5 EM  algorithms}}

{\bf \color{blue} 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 \color{dgreen} complete}' data $\bX^\tau = (\bY^\tau, \bZ^\tau)$ contain a
`{\sl \color{dgreen} missing}' component $\bZ$. The likelihood based on the
 complete data is
\[
L(\btheta; \bX) = f_\bX (\bX, \btheta).
\]

\bigskip

{\color{dblue}{\bf EM} ({\sl E}xpectation and {\sl M}aximisation) {\bf algorithms}
 }
\vspace{-8mm}
\begin{quote}
$\bullet$ {\sl \color{dred} E-step}: compute the conditional expectation
\[
Q(\btheta) = Q(\btheta | \bY, \btheta_0) \equiv E \{ \log L(\btheta; \bX) | \bY, 
\btheta_0\}
\]
$\bullet$ {\sl \color{dred} M-step}: maximise $Q(\btheta)$ to give an updated value
$\btheta_1$
\end{quote}
\vspace{-8mm}
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)$.


\newpage

{\bf 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 {\bf \color{red} 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 {\bf \color{red} 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 \underline{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$.
}

\newpage
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 {\bf \color{blue} 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).



\newpage
{\bf \color{blue} 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}
}



\newpage

{\bf Example 12}. {\color{red} 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.


\newpage

{\bf \color{dblue} 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
\vspace{-8mm}
{\color{blue}
\begin{quote}
(i) exact estimation error $\wh \btheta  - \btheta_o$ is unknown, and\\
(ii) $\wh \btheta$ is a random variable
\end{quote}}
\vspace{-8mm}
we have to gauge the error
\vspace{-8mm}
{\color{blue}
\begin{quote}
(i)
in terms of a probability average, and\\
(ii)
for all possible values of $\btheta_o \in \Theta$.
\end{quote}}
\vspace{-8mm}
Let $P_\btheta$, $E_\btheta$ and $\var_\btheta$ denote the probability distribution,
expectation and variance under $\btheta_o = \btheta$.


{\bf \color{blue} Bias}: Bias$_\btheta( \wh \btheta)= E_\btheta(\wh \btheta) - \btheta$

{\bf \color{blue} Variance}: $ \var_\btheta(\wh \btheta) $

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

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

{\bf \color{blue} 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
{\color{dgreen}
\[
\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 {\bf \color{dgreen} 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
{\color{dred} technical tractability} while the MAE leads to
estimators which is {\color{dred} more  robust} against outliers in
observations.
 



\newpage

{\bf \color{dred} 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.

\underline{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*}

\newpage

\underline{Variance of $s(\btheta)$} --- {\bf \color{red} 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 {\color{dgreen} \bf 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

{\bf Theorem 2}. {\em \color{blue}(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 {\bf a lower bound} for any
{\sl unbiased estimator} for the parameter $g(\theta)$. 
When the equality holds, $T$  
is the {\bf \color{dmagenta} minimum variance unbiased estimator (MVUE)} 
of $g(\theta)$.

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


\newpage

{\bf 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$.

\newpage
{\bf \color{dred} 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) \underline{\color{dred} 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}.  

\newpage

(b) \underline{\color{dred} 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 {\bf approximately} that
\[
\wh \btheta \sim N \big( \btheta,\; \{\calI_{X_1}(\btheta)\}^{-1}/n \big).
\]
Therefore {\sl \color{blue} asymptotically} the MLE is unbiased and attains 
the Cram\'er-Rao lower bound. Any estimator fulfilling this condition 
is called {\bf \color{dgreen} efficient}.

{\bf \color{dred} 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}$.



\newpage

{\bf \fcolorbox{blue}{yellow} {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 {\bf basic idea of bootstrap} is to adopt the above sampling procedure in
the so-called {\sl \color{red} bootstrap world}: now the population is
$f(\cdot, \wh \theta)$ which is 
\underline{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 {\bf \color{blue}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}
{\color{blue}

\psset{framearc=.4, linecolor=gray}

\begin{psmatrix}[rowsep=30pt,colsep=68pt]
Real World & {\color{dred} Bootstrap World}\\
\psframebox{
$
\begin{array}{c}
{\rm Population}\\
f(\cdot, \theta_o)\\
{\rm unknown}
\end{array} 
$} &
\psframebox{
\color{dred}
$
\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{
\color{dred}
$
\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{
\color{dred}
$
\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{\color{black}Estimating}}
\end{psmatrix} 
}
\end{center}


{\color{dgreen} 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
\vspace{-4mm}
\begin{itemize}
\item
Hall, P. (1992). The bootstrap and Edgeworth expansion. Springer, New York.
\vspace{-4mm}
\item
Efron, B. and Tibshirani, R.J. (1993). An introduction to the bootstrap. 
Chapman and Hall, New York.
\vspace{-4mm}
\item
Shao, J. and Tu, D. (1995). The jackknife and bootstrap. Springer, New York.
\vspace{-4mm}
\item
Davison, A.C. and Hinkley, D.V. (1997). Bootstrap methods and their application.
Cambridge University Press, Cambridge.
\end{itemize}
}


  






