\documentclass[12pt,space]{article}
\usepackage{amsmath}
%\usepackage{floatfig}
%\usepackage{graphicx}
%\usepackage{eepic}

%\pagestyle{empty}
\setlength{\textwidth}{17cm}
\setlength{\textheight}{9in}
\setlength{\topmargin}{-.5in}
\setlength{\oddsidemargin}{-0.5cm}

\newcommand{\A}{{\bf A}}
\newcommand{\I}{{\bf I}}
\newcommand{\bSigma}{{\bf \Sigma}}
\newcommand{\bLambda}{{\bf \Lambda}}
\newcommand{\bPhi}{{\bf \Phi}}
\newcommand{\tr}{{\rm tr}}
\newcommand{\hx}{\mbox{$\hat{x}$}}
\newcommand{\hy}{\mbox{$\hat{y}$}}
\newcommand{\bx}{{\bf x}}
\newcommand{\xbad}{x_{41}}
\newcommand{\by}{{\bf y}}
\newcommand{\bs}{{\bf s}}
\newcommand{\btheta}{{\bf \theta}}
\newcommand{\hbx}{\mbox{$\hat{\bf x}$}}
\newcommand{\hby}{\mbox{$\hat{\bf y}$}}
\newcommand{\ben}{\begin{enumerate}}
\newcommand{\een}{\end{enumerate}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\arraycolsep .1in

\begin{document}


\begin{center}
\large \bf The use of Expectation Maximization (EM) for finding the
parameters of a 2D Gaussian given missing data\\(Example 2, from DHS)
\end{center}

\vspace{.25in}

Extra notes for MAS622J/1.126J by Rosalind W. Picard

\vspace{.3in}


Let $D$ be a set of data, here containing four samples.  One sample has an element that
is either missing, known to be corrupted, or is otherwise ``bad", $D_b = x_{41} = ``*"$.
\[
D = \{ \bx_1, \bx_2, \bx_3, \bx_4 \} = \{ \left[ \begin{array}{c} 0 \\ 2 \end{array} \right], \left[ \begin{array}{c} 1 \\ 0 \end{array}\right] , \left[ \begin{array}{c} 2 \\ 2 \end{array} \right], \left[ \begin{array}{c} * \\ 4 \end{array} \right]\}
\]
A 2-D Gaussian has four parameters, $\btheta$, to estimate, for which we make the reasonable starting guess $\btheta^0$:
\[ \btheta = \left[ \begin{array}{c} \mu_1 \\ \mu_2 \\ \sigma^2_1 \\ \sigma^2_2 \end{array} \right] \; \hspace{.5in}
\btheta^0 = \left[ \begin{array}{c} 0 \\ 0 \\ 1 \\ 1 \end{array} \right]
\]
Our goal is to iteratively estimate $\btheta^i$ until its values converge.

The estimation process consists of two steps: (1) E-step and (2) M-step. In the E-step we
wish to formulate the log likelihood for all the data, marginalizing over the possible values for the unknown data. We assume there is a known set of parameters: for this, we use our current best guess of the parameter vector.

\begin{eqnarray*}
Q(\btheta; \btheta^0) & = & E_{\xbad} [ \ln p(x_g, x_b; \btheta) | D_g; \btheta^0 ] \\
 & = & \int_{-\infty}^{\infty} \left[\sum_{k=1}^3 \ln p(\bx_k | \btheta) + \ln p(\bx_4 | \btheta) \right] p(\xbad| x_{42}=4; \btheta^0)d\xbad 
% & = & \sum_{k=1}^{3}\ln p({\bf x_k}|\btheta) + \int_{- \infty}^{\infty} \ln p \left( \left[ \begin{array}{c} \xbad \\ 4 \end{array} \right] | \btheta \right) \frac{p(\xbad, x_{42} | \btheta^0)}{c} d\xbad
 \end{eqnarray*}

Note that $x_{41}$ is independent of all the other $x_{ij} \in D_g$ except for possibly $x_{42}$.

Now let $c = \int p(x_{41}, x_{42} | \btheta^0) dx_{41} = p(x_{42} |\btheta^0)$. Applying Bayes rule, we get:

\begin{eqnarray*}
Q(\btheta; \btheta^0) & = & \int_{-\infty}^{\infty} \left[\sum_{k=1}^3 \ln p(\bx_k | \btheta) + \ln p(\bx_4 | \btheta) \right] \frac{p(\xbad, x_{42}; \btheta^0)}{c}d\xbad\\ 
 & = & \sum_{k=1}^{3}\ln p({\bf x_k}|\btheta) + \int_{- \infty}^{\infty} \ln p \left( \left[ \begin{array}{c} \xbad \\ 4 \end{array} \right] | \btheta \right) \frac{p(\xbad, x_{42} | \btheta^0)}{c} d\xbad
\end{eqnarray*}




Also, since $c$ is not a function of the bad data, $x_{41}$, we can pull it out of the integral:

%Note that $c$ arises from applying Bayes rule:
%$ c = \int p(\xbad, x_{42} | \btheta^0) d\xbad  = p(x_{42} | \btheta^0)$,
%and also $c$ is not a function of the bad data, so we can pull it out of the integral:


\begin{eqnarray*}
\lefteqn{Q(\btheta; \btheta^0)}\\
&=&  \sum_{k=1}^{3}\ln p({\bf x_k}|\btheta) + \frac{1}{c} \int_{- \infty}^{\infty} \ln p \left( \left[ \begin{array}{c} \xbad \\ 4 \end{array} \right] | \btheta \right)  \frac{1}{2 \pi^{d/2} \left| \begin{array}{ll}  1 & 0 \; \\ 0 & 1 \end{array}\right| ^{1/2}} e^{-\frac{(\xbad - 0)^2}{2}} e^{-\frac{(4-0)^2}{2}}d\xbad \\
&=& \sum_{k=1}^{3}\ln p({\bf x_k}|\btheta) + \frac{1}{c} \int_{- \infty}^{\infty} \ln \left[ \frac{1}{2 \pi \left|\begin{array}{ll}  \sigma_1^2 & 0  \\ 0 & \sigma_2^2 \end{array}\right| ^{1/2}} e^{-\frac{1}{2}\left( \begin{array}{c} \xbad - \mu_1 \\ 4 - \mu_2 \end{array} \right)^T \left[ \begin{array}{cc}\sigma_1^{-2} & 0 \\ 0 & \sigma_2^{-2}\end{array} \right] \left(\begin{array}{c} \xbad - \mu_1 \\ 4 - \mu_2\end{array} \right)} \right]\\ 
& & * \left[ \frac{1}{2\pi} e^{-\frac{(\xbad^2 + 16)}{2}} \right] d\xbad
  \end{eqnarray*}

Rewrite the log term in the integral above:

\[
-\ln 2\pi\sigma_1\sigma_2 - \frac{(\xbad - \mu_1)^2}{2 \sigma_1^2} - \frac{(4-\mu_2)^2}{2\sigma_2^2}
\]

or

\[
-\ln 2\pi\sigma_1\sigma_2 - \frac{(\xbad^2 - 2\xbad\mu_1+\mu_1^2)}{2 \sigma_1^2} - \frac{(4-\mu_2)^2}{2\sigma_2^2}
\]

Now we can separate out the terms that involve the bad data, $\xbad$, and rewrite $Q$ as:

\begin{eqnarray*}
\lefteqn{Q(\btheta; \btheta^0) = }\\
&&  \sum_{k=1}^{3}\ln p({\bf x_k}|\btheta) + \left[ -\ln 2\pi\sigma_1\sigma_2 - \frac{\mu_1^2}{2\sigma_1^2} - \frac{(4-\mu_2)^2}{2\sigma_2^2}\right] \left[ \frac{1}{c}\int_{-\infty}^\infty \frac{1}{2\pi}e^{\frac{-(\xbad^2 + 16)}{2}} d\xbad \right] + III
  \end{eqnarray*}
where $III$ contains the remaining two terms involving $ -\xbad^2$ and $ 2 \xbad \mu_1$. Note first that $\frac{1}{c}$ multiplying the integral causes that whole term to go to 1.  
Subsequently,

\begin{eqnarray*}
Q(\btheta; \btheta^0) & = & I + II + III \\
& = &  \sum_{k=1}^{3}\ln p({\bf x_k}|\btheta) + \left[ -\ln 2\pi\sigma_1\sigma_2 - \frac{\mu_1^2}{2\sigma_1^2} - \frac{(4-\mu_2)^2}{2\sigma_2^2}\right] + III
\end{eqnarray*}
  
Now, let's look at that last term:

\begin{eqnarray*}
III = \frac{1}{c}\int_{-\infty}^\infty -\frac{\xbad^2}{2\sigma_1^2} \frac{1}{2 \pi}  e^{\frac{-(\xbad^2 + 16)}{2}} d\xbad + \frac{1}{c}\int_{-\infty}^\infty -\frac{2\xbad \mu_1}{2\sigma_1^2} \frac{1}{2 \pi}  e^{\frac{-(\xbad^2 + 16)}{2}} d\xbad
\end{eqnarray*}

The second integral above consists of an odd multiplier, $\frac{2 x_{41} \mu_1}{2 \sigma_1^2}$, times an even exponential (even with respect to $x_{41}$). Thus, the integral is odd and integrates to zero. This leaves only the first integral:

\begin{eqnarray*}
III &=& \frac{1}{c}\int_{-\infty}^\infty -\frac{\xbad^2}{2\sigma_1^2} \frac{1}{2 \pi}  e^{\frac{-(\xbad^2 + 16)}{2}} d\xbad\\
&=& \frac{-\frac{e^{-8}}{2\sigma_1^2}\int_{-\infty}^\infty -\frac{\xbad^2}{2 \pi} e^{-\frac{\xbad^2}{2}} d\xbad}{e^{-8} \int_{-\infty}^\infty \frac{1}{2 \pi} e^{-\frac{\xbad^2}{2}} d\xbad}
\end{eqnarray*}

Let $z = \frac{x_{41}}{\sqrt{2}}$ so $dx_{41} = \sqrt{2} d z$, and cancel the $\frac{e^{-8}}{2 \pi}$ in the numerator and denominator. Then:

\begin{eqnarray*}
III &=& \frac{-\frac{1}{2\sigma_1^2}\int_{-\infty}^\infty 2z^2 e^{-z^2} \sqrt{2} dz}{\int_{-\infty}^\infty e^{-z^2} \sqrt{2} dz}
\end{eqnarray*}

Integral tables come in handy at this point. From standard tables, we find these:

\[ \int_0^\infty x^2 e^{-x^2} dx  =  \sqrt{\frac{\pi}{16}} \; 
\hspace{.5in}
\int_0^\infty e^{-x^2} dx  =  \sqrt{\frac{\pi}{4}} 
\]



%\begin{eqnarray*}
%\int_0^\infty x^2 e^{-x^2} dx & = & \sqrt{\frac{pi}{16}}\\
%\int_0^\infty e^{-x^2} dx & = & \sqrt{\frac{\pi}{4}}
%\end{eqnarray*}


Using these relationships, we simplify:

\begin{eqnarray*}
III & = & \frac{\frac{-1}{\sigma_1^2}(2 \sqrt{2} \sqrt{\frac{\pi}{16}})}{2 \sqrt{2} \sqrt{\frac{\pi}{4}}}\\
& = & \frac{-1}{\sigma_1^2} \sqrt{\frac{4}{16}} = -\frac{1}{2 \sigma_1^2}
\end{eqnarray*}



Now, return to (1) substituting in this simplified $III$ to write Q:

\begin{eqnarray*}
Q(\btheta; \btheta^0) =  \sum_{k=1}^{3}\ln p({\bf x_k}|\btheta) -\ln 2\pi\sigma_1\sigma_2 - \frac{1 + \mu_1^2}{2\sigma_1^2} - \frac{(4-\mu_2)^2}{2\sigma_2^2}
\end{eqnarray*}

Equation (2) is called the E-step, the first step of Expectation-Maximization.

Let's expand the first term in (2) and then we will be ready for the M-step, in which we'll maximize Q:

\begin{eqnarray*}
Q(\btheta; \btheta^0) =  \sum_{k=1}^{3} \left(  -\ln 2 \pi \sigma_1 \sigma_2 - \frac{(x_{k1} -\mu_1)^2}{2\sigma_1^2} - \frac{(x_{k2}-\mu_2)^2}{2\sigma_2^2} \right)  -\ln 2\pi\sigma_1\sigma_2 - \frac{1 + \mu_1^2}{2\sigma_1^2} - \frac{(4-\mu_2)^2}{2\sigma_2^2}
\end{eqnarray*}


Now we're ready for the Maximization (M-step):
\begin{eqnarray*}
\btheta^1 = \arg\max_\btheta Q(\btheta; \btheta^0)
\end{eqnarray*}

We need to take derivatives of Q w.r.t. each of its four parameters, set the derivatives equal to zero, and solve for the new parameter values.


\begin{eqnarray*}
\frac{dQ(\btheta;\btheta^0)}{d\mu_1} & = & \frac{2}{2 \sigma_1^2} \sum_{k = 1}^3 (x_{k1} - \mu_1) - \frac{2 \mu_1}{2 \sigma_1^2} = 0\\
\sum_{k = 1}^3 x_{k1} & = &  4 \mu_1 \\
\mu_1 & = &  \frac{0 + 1 + 2}{4} = \frac{3}{4}
\end{eqnarray*}
  
  

\begin{eqnarray*}
\frac{dQ(\btheta;\btheta^0)}{d\mu_2} & = & \frac{1}{\sigma_2^2} \sum_{k = 1}^3 (x_{k2} - \mu_2) - \frac{(4- \mu_2)}{\sigma_2^2} = 0\\
\mu_2 & = &  \frac{1}{4} (2 + 0 + 2 + 4) = 2
\end{eqnarray*}
  

 

\begin{eqnarray*}
\frac{dQ(\btheta;\btheta^0)}{d\sigma_1} & = & \frac{-4(2 \pi \sigma_2)}{2 \pi \sigma_1 \sigma_2} + \sum_{k = 1}^3 \frac{(x_{k1} - \mu_1)^2}{\sigma_1^3} + \frac{(1 + \mu_1^2)}{\sigma_1^3} = 0\\
\frac{4 \sigma_1^3}{\sigma_1} & = &  (0 - \mu_1)^2 + (1 - \mu_1)^2 + (2 - \mu_1)^2 + (1 + \mu_1^2)\\
\sigma_1^2 & = & \frac{1}{4} \left(\frac{9}{16} + \frac{1}{16} +\frac{25}{16} + \frac{16}{16} +\frac{9}{16} \right) = \frac{15}{16} = 0.938
\end{eqnarray*}


Similarly, solving $\frac{dQ(\btheta;\btheta^0)}{d\sigma_2} = 0$ leads to the solution $\sigma_2^2 = 2$.

Now we have completed the M-Step and we have the next iteration of our estimate of the
parameters:

\[ \btheta = \left[ \begin{array}{c} \mu_1 \\ \mu_2 \\ \sigma^2_1 \\ \sigma^2_2 \end{array} \right] \; \hspace{.5in}
\btheta^0 = \left[ \begin{array}{c} 0 \\ 0 \\ 1 \\ 1 \end{array} \right] 
\hspace{.5in}
\btheta^1 = \left[ \begin{array}{c} 3/4 \\ 2 \\ 0.938 \\ 2 \end{array} \right]
\]

We use this to form $Q(\btheta;\btheta^1)$ and go back to the beginning of the E-step. Repeat this process until the parameters converge.

The E-M algorithm guarantees that the log-likelihood of the good data (with the bad data
marginalized) will increase monotonically.

Fun facts: E-M is a Maximum Likelihood method, not a fully Bayesian method. Also, the
Baum-Welch (forward-backward) algorithm used for training HMM's is an example of the E-M
method.

\end{document}
