\documentclass[11pt,twocolumn]{article}
\usepackage{times,amsmath,amsfonts}
\usepackage[dvips]{graphicx,graphics}

\setlength{\textheight}{10in} \setlength{\textwidth}{7.1in}
\setlength{\columnsep}{0.3in} \setlength{\topmargin}{0in}
\setlength{\headheight}{-0.3in} \setlength{\headsep}{0in}
\setlength{\oddsidemargin}{-0.6in}  % Centers text.

% Add the period after section numbers.  Adjust spacing.
%\newcommand{\Section}[1]{\vspace{-5pt}\section{\hskip -1em.~~#1}\vspace{-3pt}}
%\newcommand{\SubSection}[1]{\vspace{-2pt}\subsection{\hskip -1em.~~#1}
%        \vspace{-3pt}}
\newcommand{\bqn}{\begin{eqnarray*}}
\newcommand{\eqn}{\end{eqnarray*}}
\newcommand{\bq}{\begin{eqnarray}}
\newcommand{\eq}{\end{eqnarray}}

\begin{document}

\title{Stat 992: Lecture 32\\
Gaussian mixture model I.}
\author{Moo K. Chung \tt{mchung@stat.wisc.edu}}
\maketitle \thispagestyle{empty}

%\begin{figure}
%\centering
%\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.8]{homework05-1.eps}
%\end{figure}

\begin{enumerate}


\item {\em Gaussian mixture model.} Download data and run the
following {\tt MATLAB} code.
\begin{verbatim}
load sagittal.data;
intensity=reshape(sagittal,256*124,1);
intensity=intensity(find(
(intensity >=50)&(intensity<=250)));
hist(intensity,100);
\end{verbatim}

We will get 19,660 values for $\tt{intensity}$ which are gray
scale values between 50 and 250 of the magnetic resonance image
(MRI) of the human brain. Let $f(y_i)$ be the image intensity at
pixel position $y_i$. From the intensity histogram, we see that
the image intensity can be modelled as a two component Gaussian
mixture
$$f(y) = p f_1(y) + (1-p)f_2(y)$$
where $f_1 \sim N(\mu_1,\sigma_1^2)$ and $f_2 \sim
N(\mu_2,\sigma_2^2)$. We will estimate parameters by maximizing
the likelihood function using the EM algorithm.

\item {\em Expectation-maximization} (EM) method is an iterative
method for maximizing  difficult likelihood problems. It was first
introduced by Dempster {\em et al.} (J. Roy. Statist. Soc. 1997).
Read Robert and Casella's Monte Carlo Statistical Methods for the
introduction to EM. Flury's A First Course in Multivariate
Statistics for the detailed discussion on EM applied to Gaussian
mixture model. See Little and Rubin (1987) and McLachlan and
Krishnan (1997). EM is widely used in imaging. Instead of
classifying pixels into classes, it is possible to assign the
probability belong to classes (SPM'99 segmentation algorithm).

Suppose we have a random sample $X=\{X_1,\cdots, X_n\} \sim
f(x|\theta)$. We wish to find the maximum likelihood estimator of
$\theta$:
$$\hat \theta = \arg \max_{\theta} \prod_{i=1}^n f(x_i|\theta) = \arg \max_{\theta} \sum_{i=1}^n \ln f(x_i|\theta).$$
This optimization problem is nontrivial when $f$ is complex and it
gives a motivation for the development of the EM algorithm. We
augment the data with with {\em latent} (unobserved) data $X^m$
such that the complete data $X^c = (X, X^m)$. The density of $X^c$
is denoted as
$$X^c=(X,X^m) \sim f(x^c) = f(x,x^m).$$ The conditional density for the missing data $X^m$ is
$$f(x^m|x,\theta) = \frac{f(x,x^m|\theta)}{f(x|\theta)}.$$
Rearranging terms, $f(x|\theta) =
\frac{f(x,x^m|\theta)}{f(x^m|x,\theta)}.$ Taking logarithm, we get
log-likelihood for the observed data $$\ln f(X|\theta) = \ln
f(X^c|\theta) - \ln f(X^m|X,\theta).$$ Take expectation with
respect to $f(x^m|x,\theta_0)$ on both sides so that the left side
is taken as constant.
\bqn \ln f(X|\theta) &=& \mathbb{E} [\ln f(X^c|\theta) |X,\theta_0]\\
&& - \mathbb{E} [\ln f(X^m|X,\theta)|X,\theta_0].\eqn


\begin{figure}
\centering
%\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.6]{homework05-1.eps}
\caption{Left: mid sagittal brain image. Right: histogram of image
intensity.}
\end{figure}

\item Following Robert and Casella, denote the expected
log-likelihood for the complete data by
$$Q(\theta|\theta_0,X) = \mathbb{E} [\ln f(X^c|\theta)|X,\theta_0].$$
We maximize the likelihood in the following fashion.\\
1. E-step: compute $Q(\theta|\hat \theta_{j-1},X)$.\\
2. M-step: maximize $Q(\theta|\hat \theta_{j-1},X)$ and take
$$\hat \theta_j = \arg \max_{\theta} Q(\theta|\hat \theta_{j-1},X).$$
If the above procedure is iterated, we get the sequence of
estimators $\hat \theta_0,\hat \theta_1, \cdots $ and it can be
shown that it converges to
the maximum likelihood estimator $\hat \theta$.\\
{\em Proof.} By the definition of $\hat \theta_{j+1}$, $$Q(\hat
\theta_{j+1}|\hat \theta_j,x) \geq Q(\hat \theta_j|\hat
\theta_j,x).$$ Let $R(\theta|\theta_0,X) = \mathbb{E} [\ln
f(X^m|X,\theta)|X,\theta_0].$ It can be shown that \bq\theta_0 =
\arg \max_{\theta} R(\theta|\theta_0,x) \label{eq:jensen}\eq

Then $ R(\hat \theta_{j+1}|\hat \theta_j,x) \leq R(\hat
\theta_j|\hat \theta_j,x).$ Consequently \bqn\ln f(x|\hat
\theta_j) &=& Q(\hat \theta_j|\hat \theta_j,x) - R(\hat
\theta_j|\hat\theta_j,x)\\
&\leq& Q(\hat \theta_{j+1}|\hat \theta_j,x) - R(\hat
\theta_{j+1}|\hat\theta_j,x)\\
 &\leq& \ln f(x|\hat \theta_{j+1}).\eqn This inequality
guarantees the the sequence of estimators $\hat \theta_j$
monotonically increase the likelihood. Since the monotonically
increasing sequence is bounded, i.e. $\ln f(x|\hat \theta_j) \leq
\ln f(x|\hat \theta)$, it must be converging to a stationary point
but it is not clear if the limit is $\ln f(x|\hat \theta)$. To
guarantee that the limit converges to the maximum likelihood
estimator, additional conditions are needed. See Boyles (J. Roy.
Statisti. Soc. B. 1983) and Wu (Ann. Statist. 1983)

The difficulty of EM algorithm is at the E-step where we need to
compute the conditional expectation $Q(\theta|\hat
\theta_{j-1},x)$. The {\em Monte-Carlo EM} algorithm overcome this
by simulating missing data $X^m \sim f(x^m|x,\theta)$. so that
$$\widehat Q(\theta|\theta_0,x) = \frac{1}{l}\sum_{j=1}^l \ln
f(X,X^m|\theta).$$ {\bf Problem 46.} Prove inequality
(\ref{eq:jensen}). Hint. Use Jensen's inequality. \item {\em
Multinomial example.} $X =(X_1, \cdots, X_n) \sim
M(n,p_1,p_2,\cdots, p_{k-1},p_k)$ is multinomial if
$$P(X_1=x_1, \cdots, X_k=x_k) = \frac{n!}{x_1 ! \cdots x_k !}p_1^{x_1} \cdots p_k^{x_k},$$
where $\sum_{i=1}^k x_i = n$ and $\sum_{i=1}^k p_i =1$. The
following example is motivated from Rao C.R. (Linear statistical
inference and its applications, 1973) and simplified from $k=4$ to
$k=3$. Suppose observations $x=(x_1,x_2,x_3)$ are from the
multinomial $M(n,\frac{1}{2}+\frac{p}{4},\frac{1-p}{2}
,\frac{p}{4})$. Find the MLE of $p$.\\
{\em Solution.} The observed likelihood is
$$L(p|x) \propto (2+p)^{x_1}(1-p)^{x_2}p^{x_3}.$$ MLE can be computed if
 one is willing to do algebraic manipulation. Instead we augment data by additional data such that
$$x^c =(z_1,z_2,x_2,x_3) \sim M(n,\frac{1}{2},\frac{p}{4},\frac{1-p}{2} ,\frac{p}{4})$$
with $x_1=z_1+z_2$. The complete likelihood function is $L(p|x^c)
\propto p^{z_2+x_3}(1-p)^{x_2}.$ This is easier to manipulate than
$L(p |x)$. The conditional expected complete log-likelihood $Q(p
|p_0,x)$ is
$$ \mathbb{E}[\log L(p|X^c)|p_0,x] \propto \mathbb{E}_{p_0} [(Z_2+x_3)\log p + x_2 \log(1-p)]
$$
We can also show that $Z_2|X_1=x_1 \sim B \Big(x_1, \frac{p}{2+p}
\Big)$ (check it yourself). So the E-step is
 $$Q(p|p_0,x) = \Big(\frac{p_0}{2+p_0}x_1+x_3\Big)\log p + x_2 \log(1-p).$$
The M-step is also easy. Computing $\partial Q/\partial p =0$, we
get
$$\hat p_1 =
\frac{\frac{p_0}{2+p_0}x_1+x_3}{\frac{p_0}{2+p_0}x_1+x_2+x_3}.$$
This gives a recursive formula to find $\hat p_i$.
\begin{verbatim}
clear p; x1=35;x2=50;x3=16;
p(1)=0.9;
for i=2:20
  pr= p(i-1)/(2+p(i-1))
  p(i) = (pr*x1 + x3)/(pr*x1+x2+x3);
end; hold on; plot(p)
\end{verbatim}
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.24]{lecture27-1.eps}
\includegraphics[scale=0.24]{lecture27-2.eps}
\caption{Left: convergence of EM algorithm. Right: Monte-Carlo EM
algorithm.}
\end{figure}
\end{enumerate}
\end{document}

\item Consider following random field model on manifolds $Y_i(x) =
\mu_i(x) +\epsilon(x), x \in \partial \Omega$ for subject $j$. We
assume i.i.d. zero mean Gaussian random field $\epsilon_j$ with
isotropic covariance $R(x,y) = \sigma^2f(d(x,y))$ where $d$ is a
geodesic distance and $f$ some function.

Suppose there are other measurement $Z_j \sim N(\mu_Z
,\sigma_Z^2)$. The sample correlation between two measurements
$Y_j$ and $Z_j$ at fixed point $x$ is given by
$$\rho(x) = \frac{\mbox{Cov}(Y(x),Z)}{$$


\item {\em Covariance function estimation.} On manifolds $\partial
\Omega$, we have model
$$Y(x) = \mu(x) + \epsilon(x)$$
where $\epsilon$ is unknown random field. We may estimate $\mu$ by
kernel smoothing so that $\widehat{\mu} = K_{\sigma}*Y$. Now we
estimate the noise with the residual $\widehat{\epsilon} = Y -
K_{\sigma}*Y$.