\documentclass[10pt,final]{article}
\usepackage{graphics,amsmath}
\usepackage{graphicx}
\usepackage{algorithmic}
\usepackage[boxed]{algorithm}
\usepackage{subfigure,mydefs}

%Comment out this line when document is finalized
%\usepackage[outline,bottom,first]{draftcopy}


\setlength{\parindent}{0cm}
\setlength{\parskip}{1.2ex plus0.2ex minus0.1ex}

\addtolength{\topmargin}{-0.5in}
\addtolength{\textheight}{1.0in}

\addtolength{\oddsidemargin}{-0.4in}
\addtolength{\evensidemargin}{-0.4in}
\addtolength{\textwidth}{0.9in}
\addtolength{\marginparwidth}{-0.4in}

\renewcommand{\thefootnote}{\fnsymbol{footnote}}

\title{\textbf{Using EM To Estimate A Probability Density With A Mixture Of Gaussians}}
\author{\textbf{Aaron A. D'Souza}\\\small\textsf{adsouza@usc.edu}}
%\\\small\textsf{http://www-slab.usc.edu/\~{}adsouza}
\date{}

% \newcommand{\Zim}{z_{ik}} 
\newcommand{\Zim}{\pi_{ik}}
\newcommand{\mSigmam}{\mSigma_k}
\newcommand{\imSigma}{\mSigma^{-1}}
\newcommand{\imSigmam}{\mSigma^{-1}_k}
\newcommand{\pim}{\pi_k}

\allowdisplaybreaks[3]

\begin{document}
\maketitle

\section{Introduction}
\label{sec:introduction}


The problem we are trying to address in this note is simple. Given a
set of data points $\mX = \Brace{\vx_i}^N_1$, we wish to determine the
underlying probability distribution $p(\vx)$, that generates this data.
In general, the distribution $p(\vx)$ could be any real-valued, scalar
function with the following constraints:
\[
p(\vx) \geq 0, \forall \vx \text{ and } \int_{-\infty}^{\infty} \!\!\!\!\!p(\vx) d\vx = 1
\]

The framework of mixture modeling using Gaussians makes the following
assumptions:

\begin{itemize}
\item The data was generated using a \emph{set} of $K$ probability
  distributions.
  
\item Each of the individual probability distributions is a Gaussian:
\[
 \vx;\vtheta_k \sim \Normdist{\vx}{\vmu_k}{\mSigmam}\footnotemark
\]
\footnotetext{We use the notation $\Normdist{\vx}{\vmu}{\mSigma}$ to denote the
  multivariate Gaussian distribution which is mathematically defined as:
  \begin{equation*}
    \Normdist{\vx}{\vmu}{\mSigma} = (2\pi)^{-d/2}\Det{\mSigma}^{-1/2}\exp
  \Brace{-\frach(\vx-\vmu)^T\imSigma(\vx-\vmu)}
  \end{equation*}
  where $d$ is the dimensionality of $\vx$.}
i.e., the probability of generating a data point $\vx$ under the
$m^{th}$ model is given according to a Gaussian distribution with mean
$\vmu_k$ and covariance $\mSigmam$.

\item Each data point is generated according to the following
  algorithm:\\
  \begin{algorithmic}[1]
    \FOR{$i=1$ to $N$}
    \STATE $m \leftarrow$ index of one of the $K$ models randomly selected\\according to the prior probability vector $\vpi$
    \STATE Randomly generate $\vx_i$ according to the distribution $\Normdist{\vx_i}{\vmu_k}{\mSigmam}$
    \ENDFOR
  \end{algorithmic}
\end{itemize}

Representing probabilistic systems as graphical models is rapidly
becoming a useful tool in Bayesian analysis. The graphical model
corresponding to our formulation of the data generation process is
shown in fig.~\ref{fig:graphical_model}. We adopt the convention that
circular nodes correspond to random variables in the model, while
rectangular nodes correspond to variables that parameterize the
distributions of these random variables. Circular nodes with double
borders indicate observed variables in our model (the data).

% \begin{figure}[t]
%   \begin{center}
%     \scalebox{0.7}{\includegraphics{graphical_kodel.eps}}
%     \caption{Graphical model for maximum likelihood density
%       estimation using a mixture of Gaussians}
%     \label{fig:graphical_kodel}
%   \end{center}
% \end{figure}

\begin{figure}[t]
  \begin{center}
    \scalebox{0.85}{\input{figures/graphical_model.pstex_t}}
    \caption{Graphical model for maximum likelihood density
      estimation using a mixture of Gaussians}
    \label{fig:graphical_model}
  \end{center}
\end{figure}

In this model we have introduced an additional variable $\vz_i$
associated with each $\vx_i$.  The $\vz_i$ variables are \emph{indicator}
variables that are multinomially distributed according to the
parameter vector $\vpi$, and indicate which component generates the
corresponding $\vx_i$. It is easiest to think of each $\vz_i$ as an $K$
dimensional vector with a 1 in the element corresponding to the
selected mixture component, and 0's in all other elements. The
probability of the $m^{th}$ element being 1, is $w_k$.

\[
 \vz_i = \underbrace{\BMatrix{0 & 0 & \cdots & 0 & 1 & 0 & \cdots & 0 & 0}^T}_{\text{$K$ elements}}
\]

Since we do not know the corresponding $\vz_i$ for each $\vx_i$ (if we
did, then we simply group the $\vx_i$ according to their $\vz_i$, and
fit a single Gaussian to each group), these variables are called
\emph{hidden} variables.

Our problem is now reduced to finding the values of the model
parameters $\vmu_k$ and $\mSigmam$ for each of the $K$ models, as well
as the prior probability vector $\vpi$, which when plugged into the
generative model, is most likely to generate the observed data
distribution. In other words we are interested in maximizing the
\emph{likelihood} $\mathcal{L}(\vtheta)=p(\mX;\vtheta)$ of generating
the observed data given the model parameters $\vtheta =
\Brace{\vmu_k,\mSigmam,w_k}_1^K$. This approach is called the
Maximum Likelihood (ML) framework since it finds the parameter
settings that maximize the likelihood of observing the data.

Although the ML approach is an intuitively appealing solution, we
often find that maximizing the expressions for likelihood w.r.t. the
parameters $\vtheta$ are often analytically intractable. The
Expectation Maximization (EM) algorithm can be used to simplify the
math considerably.

\begin{figure}[t]
  \begin{center}
  \includegraphics[width=0.8\textwidth,natwidth=610,natheight=742]{mixtures_drawing.ps}
    \caption{Mixtures equations}
    \label{fig:mixtures}
  \end{center}
\end{figure}


\newpage
\section{Estimating the Model Parameters using EM}
\label{sec:em}

Instead of attempting to maximize the likelihood of the observed data
$p(\mX;\vtheta)$, we attempt instead to maximize the likelihood of the
joint distribution of $\mX$ and $\mZ = \Brace{\vz_i}_1^N$,
$p(\mX,\mZ;\vtheta)$. For the purposes of maximization we can also
work with the logarithm of this quantity, $\log p(\mX,\mZ;\vtheta)$.
This quantity is also known as the \emph{complete} log-likelihood.
Since we cannot observe the values of the random variables $\vz_i$ we
must work with the expectation of this quantity w.r.t. some
distribution $Q(\mZ)$.
%  From Eq.~\eqref{eq:lower_bound} we can see that
% for any arbitrary distribution $Q(\mZ)$, maximizing the expectation of
% the complete log-likelihood $\Expect{\log
%   p(\mX,\mZ|\vtheta)}_{Q(\mZ)}$, is equivalent to maximizing a lower
% bound to the incomplete data likelihood. This in itself is useless
% since we are only maximizing a lower bound and have no guarantee of
% maximizing $p(\mX|\vtheta)$. However, using Eq.~\eqref{eq:jensen} we
% can also show that for the special case where
% $Q(\mZ)=p(\mZ|\mX,\vtheta)$, i.e. when the expectation of the complete
% log-likelihood is taken w.r.t. the posterior distribution of the
% hidden variables, the lower bound becomes an equality, and we are
% guaranteed to maximize the incomplete log-likelihood $p(\mX|\vtheta)$.

The log of the complete data likelihood can be written as follows:
\begin{align*}
  l_c(\vtheta) & = \log p(\mX,\mZ;\vtheta)\\
  & = \log \prod_i^N p(\vx_i,\vz_i;\vtheta)\\
  & = \log \prod_i^N \prod_k^K \Bracket{p(\vx_i|\Zim=1;\vtheta)p(\Zim=1)}^{\Zim}\\
  & = \sumin \sum_k^K \Zim \log p(\vx_i|\Zim=1;\vtheta) + \Zim \log w_k
\end{align*}

Since we have assumed that each of the individual models is a
Gaussian, the quantity $p(\vx_i|m,\vtheta)$ is simply the conditional
probability of generating $\vx_i$ given that the $m^{th}$ model is
chosen:\\
\begin{equation}
  \label{eq:cond_prob}
  \log p(\vx_i|\Zim=1;\vtheta) = \frac{1}{\Paren{2\pi}^{d/2}\Det{\mSigmam}^{1/2}}\exp
  \Brace{-\frach \Paren{\vx_i-\vmu_k}^T\imSigmam\Paren{\vx_i-\vmu_k}}
\end{equation}

Taking expectations w.r.t. $Q(\mZ)$ we get:\\
\begin{equation}
  \label{eq:expect_log_likelihood}
  \Expect{l_c(\vtheta)}_{Q(\mZ)} = \sumin \sum_k^K \Expect{\Zim} \log p(\vx_i|\Zim=1;\vtheta) + \Expect{\Zim} \log w_k
\end{equation}


\subsection{The M step}
\label{sec:m_step}

The ``M'' step in EM takes the expected complete log-likelihood as
defined in eq. \eqref{eq:expect_log_likelihood} and maximizes it
w.r.t. the parameters that are to be estimated; in this case $w_k$,
$\vmu_k$, and $\mSigmam$.

Differentiating eq. \eqref{eq:expect_log_likelihood} w.r.t. $\vmu_k$ we get:
\begin{equation}
  \label{eq:diff_ku}
  \frac{\partial \Expect{l_c(\vtheta)}_{Q(\mZ)}}{\partial \vmu_k} = 
  \sumin \Expect{\Zim} \frac{\partial}{\partial \vmu_k}\log p(\vx_i|\Zim=1;\vtheta) = \vzero
\end{equation}

We can compute $\frac{\partial}{\partial \vmu_k}\log p(\vx_i|\Zim=1;\vtheta)$
using eq. \eqref{eq:cond_prob} as follows:
\begin{align*}
  \frac{\partial}{\partial \vmu_k}\log p(\vx_i|\Zim=1;\vtheta) & =
  \frac{\partial}{\partial \vmu_k}\log
  \Brace{\frac{1}{\Paren{2\pi}^{d/2}\Det{\mSigmam}^{1/2}}\exp
    \Brace{-\frach \Paren{\vx_i-\vmu_k}^T\imSigmam\Paren{\vx_i-\vmu_k}}}\\
  & = - \frach \frac{\partial}{\partial \vmu_k} \Paren{\vx_i-\vmu_k}^T\imSigmam\Paren{\vx_i-\vmu_k}\\
  & = \Paren{\vx_i-\vmu_k}^T\imSigmam\footnotemark
\end{align*}
\footnotetext{Where we have used the relation $\frac{\partial}{\partial \vx}\vx^T\mA\vx = \vx^T(\mA+\mA^T)$}

Substituting this result into eq. \eqref{eq:diff_ku}, we get:
\[
  \sumin\Expect{\Zim}\Paren{\vx_i-\vmu_k}^T\imSigmam = \vzero\\
\]  
giving us the update equation:
\begin{equation}
  \label{eq:mu_update}
  \vmu_k = \frac{\sumin\Expect{\Zim}\vx_i}{\sumin\Expect{\Zim}}
\end{equation}

Differentiating eq. \eqref{eq:expect_log_likelihood} w.r.t. $\imSigmam$ we get:
\begin{equation}
  \label{eq:diff_sigma}
  \frac{\partial \Expect{l_c(\vtheta)}_{Q(\mZ)}}{\partial \imSigmam} = 
  \sumin \Expect{\Zim} \frac{\partial}{\partial \imSigmam}\log p(\vx_i|\Zim=1;\vtheta) = \vzero
\end{equation}

We can compute $\frac{\partial}{\partial \imSigmam}\log p(\vx_i|\Zim=1;\vtheta)$ using eq. \eqref{eq:cond_prob} as follows:
\begin{align*}
  \frac{\partial}{\partial \imSigmam}\log p(\vx_i|\Zim=1;\vtheta) & =
  \frac{\partial}{\partial \imSigmam}\log
  \Brace{\frac{1}{\Paren{2\pi}^{d/2}\Det{\mSigmam}^{1/2}}\exp
    \Brace{-\frach \Paren{\vx_i-\vmu_k}^T\imSigmam\Paren{\vx_i-\vmu_k}}}\\
  & = \frac{\partial}{\partial \imSigmam} \Brace{\frach\log\Det{\imSigmam} - \frach\Paren{\vx_i-\vmu_k}^T\imSigmam\Paren{\vx_i-\vmu_k}}\\
  & = \frach \mSigmam - \frach\Paren{\vx_i-\vmu_k}\Paren{\vx_i-\vmu_k}^T\footnotemark
\end{align*}
\footnotetext{Where we have used the relation $\frac{\partial}{\partial \mX}\log\Det{\mX} = \Paren{\mX^{-1}}^T$ and $\frac{\partial}{\partial \mA}\vx^T\mA\vx = \vx\vx^T$}

Substituting this result into eq. \eqref{eq:diff_sigma}, we get:
\[
  \sumin\Expect{\Zim}\Paren{\frach \mSigmam - \frach\Paren{\vx_i-\vmu_k}\Paren{\vx_i-\vmu_k}^T} = \vzero
\]

giving us the update equation:
\begin{equation}
  \label{eq:sigma_update}
  \mSigmam = \frac{\sumin\Expect{\Zim}\Paren{\vx_i-\vmu_k}\Paren{\vx_i-\vmu_k}^T}{\sumin\Expect{\Zim}}
\end{equation}

In order to maximize the expected log-likelihood in eq.
\eqref{eq:expect_log_likelihood} w.r.t. $w_k$, we have to keep in
mind that the maximization has the constraint that $\sum_k^Kw_k =
1$. In order to enforce this constraint we use the Lagrange multiplier
$\lambda$, and augment eq.  \eqref{eq:expect_log_likelihood} as
follows:
\begin{equation}
  \label{eq:aug_exp_log_lik}
  L'\Paren{\vtheta} = \Expect{l_c(\vtheta)}_{Q(\mZ)} - \lambda\Paren{\sum_k^Kw_k - 1}
\end{equation}

We now differentiate this new expression w.r.t. each $w_k$ giving us:
\[
  \frac{\partial}{\partial w_k}\Expect{l_c(\vtheta)}_{Q(\mZ)} - \lambda = 0 \qquad \text{for } 1\leq m\leq K
\]

Using eq. \eqref{eq:expect_log_likelihood} we get:
\begin{equation}
  \left.
    \begin{aligned}
      \frac{1}{w_k}\sumin\Expect{\Zim} - \lambda & = 0\\
      \text{or equivalently }\sumin\Expect{\Zim} - \lambda w_k & = 0 \label{eq:diff_pi}
    \end{aligned}
  \:\right\}\qquad \text{for } 1\leq m\leq K
\end{equation}

Summing eq. \eqref{eq:diff_pi} over all K models we get:
\[
  \sum_k^K\sumin \Expect{\Zim} - \lambda \sum_k^K w_k = 0
\]

But since $\sum_k^Kw_k = 1$ we get:
\begin{equation}
  \label{eq:lambda}
  \lambda = \sum_k^K\sumin \Expect{\Zim} = N
\end{equation}

Substituting this result back into eq. \eqref{eq:diff_pi} we get the
following update equation:
\begin{equation}
  \label{eq:pi_update}
  w_k = \frac{\sumin\Expect{\Zim}}{N}
\end{equation}

which preserves the constraint that $\sum_k^Kw_k = 1$.












\subsection{The E step}
\label{sec:e_step}

Now that we have derived the update equations that maximize the
expected \emph{complete} log-likelihood $\Expect{\log
  p(\mX,\mZ;\vtheta)}$, we wish to ensure that we are indeed also
maximizing the \emph{incomplete} log-likelihood $p(\mX;\vtheta)$ (which
is the quantity that we are truly interested in maximizing).

As we mentioned earlier in section \ref{sec:em}, we are guaranteed to
maximize the incomplete log-likelihood only when the expectation is
taken w.r.t. the posterior distribution of $\mZ$, namely
$p(\mZ|\mX;\vtheta)$. Hence each of the expectations $\Expect{\Zim}$
that appear in the update equations derived in the previous section
(section \ref{sec:m_step}), should be computed as follows:
\begin{align*}
  \label{eq:expectation}
  \Expect{\Zim}_{p(\mZ|\mX;\vtheta)} & = 1\cdot p(\Zim=1|\vx_i;\vtheta) + 0\cdot p(\Zim=0|\vx_i;\vtheta)\\
  & = p(\Zim=1|\vx_i;\vtheta)\\
  & = \frac{p(\vx_i|\Zim=1;\vtheta)p(\Zim=1)}{\sum_{k=1}^K p(\vx_i|z_{ik}=1;\vtheta)p(z_{ik}=1)}\\
  & = \frac{p(\vx_i|\Zim=1;\vtheta)w_k}{\sum_{k=1}^K p(\vx_i|\Zim=1;\vtheta)w_k}\\
  & = \frac{p(\vx_i|\Zim=1;\vtheta)w_k}{notmalized\_over\_k}
\end{align*}















\section{Summary}
\label{sec:summary}

\renewcommand{\algorithmiccomment}[1]{\textit{//#1}}
\renewcommand{\algorithmicrequire}{\textbf{Initialize:}}
\begin{algorithm}[t]
\caption{Estimate $w_k, \vmu_k, \mSigma_k$ for $1\leq m\leq K$}
\label{algo:algorithm}
\begin{algorithmic}[1]
  \REQUIRE all $\Expect{\Zim}$, $w_k$, $\vmu_k$, and $\mSigmam$
  \REPEAT
  \FOR[The E step]{$i=1$ to $N$}
  \FOR{$m=1$ to $K$}
  \STATE \begin{align*}
    p(\vx_i|\Zim=1;\vtheta) & = (2\pi)^{-d/2}\Det{\mSigmam}^{-1/2}\exp\Brace{-\frach(\vx_i-\vmu_k)^T\mSigmam^{-1}(\vx_i-\vmu_k)}\\
    \Expect{\Zim} & = \frac{p(\vx_i|\Zim=1;\vtheta)w_k}{normalized\_over\_k}
  \end{align*}
  \ENDFOR
  \ENDFOR
  \FOR[The M step]{$m=1$ to $K$}
  \STATE \begin{align*}
    \mSigmam & = \frac{\sumin\Expect{\Zim}\Paren{\vx_i-\vmu_k}\Paren{\vx_i-\vmu_k}^T}{\sumin\Expect{\Zim}}\\
    \vmu_k & = \frac{\sumin\Expect{\Zim}\vx_i}{\sumin\Expect{\Zim}}\\
    w_k & = \frac{\sumin\Expect{\Zim}}{N}
  \end{align*}
  \ENDFOR
  \UNTIL{model parameters converge}
\end{algorithmic}
\end{algorithm}

To summarize, given a set of data points $\mX$, if we wish to estimate
the underlying probability distribution using EM to fit $K$ Gaussians,
we apply Algorithm \ref{algo:algorithm}, iterating until convergence of
the model parameters.

\section{Why does this work? (A brief review of EM theory)}
\label{sec:em_theory}

Knowing that $\log p(\mX;\vtheta)$ is difficult to maximize
analytically, we (seemingly arbitrarily) chose to maximize the
expected \emph{complete} log-likelihood $\Expect{\log
  p(\mX,\mZ;\vtheta)}_{Q(\mZ)}$ in the hope that this also increases
the \emph{incomplete} log-likelihood $\log p(\mX;\vtheta)$ (the
quantity we are \emph{really} interested in). This section will
justify this choice and prove that we are indeed maximizing $\log
p(\mX;\vtheta)$

Let us rewrite $\log p(\mX|\vtheta)$ as follows:

\begin{align}
  \log p(\mX;\vtheta) & = \log \int p(\mX,\mZ;\vtheta) d\mZ\nonumber\\
  & = \log \int Q(\mZ)\frac{p(\mX,\mZ;\vtheta)}{Q(\mZ)} d\mZ\\
  & \geq \int Q(\mZ)\log \frac{p(\mX,\mZ;\vtheta)}{Q(\mZ)} d\mZ
         \qquad \text{(Jensen's inequality)} \label{eq:jensen}\\
  & = \int Q(\mZ)\log p(\mX,\mZ;\vtheta) d\mZ - \int Q(\mZ)\log Q(\mZ) d\mZ\\
  \intertext{Hence we have the following lower-bound to $\log p(\mX;\vtheta)$:}
  \log p(\mX;\vtheta) & \geq \underbrace{\Expect{\log p(\mX,\mZ;\vtheta)}_{Q(\mZ)}}_{\text{exp. comp. log-lik.}} + \underbrace{\calH\Bracket{Q(\mZ)}}_{\text{entropy of }Q(\mZ)} = \calF(Q,\vtheta)
  \label{eq:lower_bound}
\end{align}

Since $Q(\mZ)$ is an arbitrary distribution, it is independent of
$\vtheta$. Hence in order to maximize the functional
$\calF(Q,\vtheta)$ w.r.t.~$\vtheta$, it suffices to simply maximize
$\Expect{\log p(\mX,\mZ;\vtheta)}_{Q(\mZ)}$. (Hence the M-step).

Does this maximization achieve our aim? Eq.~\eqref{eq:lower_bound}
shows that the functional $\calF(Q,\vtheta)$ is a \emph{lower bound}
to the quantity we are interested in. In which case, maximizing
$\calF(Q,\vtheta)$ does not guarantee that we are improving $\log
p(\mX;\vtheta)$ at all! If however, we set $Q(\mZ) =
p(\mZ|\mX;\vtheta)$ in Eq.~\eqref{eq:jensen}, then we see that
\emph{the lower bound in fact becomes an equality}.

\begin{align*}
  \int Q(\mZ)\log \frac{p(\mX,\mZ;\vtheta)}{Q(\mZ)} d\mZ & = \int p(\mZ|\mX;\vtheta)\log \frac{p(\mX,\mZ;\vtheta)}{p(\mZ|\mX;\vtheta)} d\mZ\\
  & = \int p(\mZ|\mX;\vtheta)\log \frac{p(\mZ|\mX;\vtheta)p(\mX;\vtheta)}{p(\mZ|\mX;\vtheta)} d\mZ\\
  & = \int p(\mZ|\mX;\vtheta)\log p(\mX;\vtheta) d\mZ\\
  & = \log p(\mX;\vtheta) \int p(\mZ|\mX;\vtheta) d\mZ\\
  & = \log p(\mX;\vtheta)
\end{align*}

This means that when computing the expected complete log-likelihood
$\Expect{\log p(\mX,\mZ;\vtheta)}_{Q(\mZ)}$, the expectation should be
taken w.r.t. the true posterior $p(\mZ|\mX;\vtheta)$ of the hidden
variables (the E step).

\section{Examples}
\label{sec:examples}

Figures \ref{fig:gauss2}, \ref{fig:gauss3} and \ref{fig:gauss4} show
the result of fitting 2, 3, and 4 Gaussians respectively to a set of
data points. In each figure the first plot shows the positions of the
means and the relative covariances of each Gaussian, while the second
shows the resulting estimated distribution obtained by marginalizing
over the models as follows:
\[
p(\vx;\vtheta) = \sum_k^K p(\vx|z_k=1;\vtheta)p(z_k=1)
\]

Figure \ref{fig:gauss5} demonstrates that one can model rather
arbitrary non-gaussian distributions provided we have a sufficient
number of mixture components. The tradeoff here is that too few
components will fail to model the structure of the data, while too
many will ``overfit'' the data. The \emph{model selection} problem in
this context is determining an appropriate compromise.

\begin{figure}[t]
  \begin{center}
    \parbox{\columnwidth}{\centering
      \subfigure{
        \scalebox{0.9}{\includegraphics{figures/mahal2.eps}}
        }\\
      \subfigure{
        \scalebox{0.9}{\includegraphics{figures/dens2.eps}}
        }
      }
    \caption{Fitting with 2 Gaussians}
    \label{fig:gauss2}
  \end{center}
\end{figure}

\begin{figure}[t]
  \begin{center}
    \parbox{\columnwidth}{\centering
      \subfigure{
        \scalebox{0.9}{\includegraphics{figures/mahal3.eps}}
        }\\
      \subfigure{
        \scalebox{0.9}{\includegraphics{figures/dens3.eps}}
        }
      }
    \caption{Fitting with 3 Gaussians}
    \label{fig:gauss3}
  \end{center}
\end{figure}

\begin{figure}[t]
  \begin{center}
    \parbox{\columnwidth}{\centering
      \subfigure{
        \scalebox{0.9}{\includegraphics{figures/mahal4.eps}}
        }\\
      \subfigure{
        \scalebox{0.9}{\includegraphics{figures/dens4.eps}}
        }
      }
    \caption{Fitting with 4 Gaussians}
    \label{fig:gauss4}
  \end{center}
\end{figure}

\begin{figure}[t]
  \begin{center}
    \parbox{\columnwidth}{\centering
      \subfigure{
        \scalebox{0.6}{\includegraphics{figures/mahal5.eps}}
        }\\
      \subfigure{
        \scalebox{0.7}{\includegraphics{figures/dens5.eps}}
        }
      }
    \caption{A sufficient number of mixture components can model arbitrary distributions}
    \label{fig:gauss5}
  \end{center}
\end{figure}

% \begin{appendix}

% \section{Sample Matlab code}
% \label{sec:matlab_code}

% The following Matlab code is a simple implementation of density
% estimation using a mixture of Gaussians.

% \scriptsize\begin{verbatim}
% function [mu, sigma, priors] = estimate_kixture_components(X, M)

% % Get dimensionality d and number of data points N
% [d,N] = size(X);

% % Compute the variance along each axis of input space
% C = diag(cov(X'));

% % Randomly initialize the means to the mean of the data
% % but randomly perturb the mean of each model by an amount
% % that's proportional to the 1/10th the std-deviation in
% % each dimension
% mu = repmat(mean(X')',1,M) + sqrt(diag(C))*randn(d,M)/10;

% % Initialize the covariance matrices to be axis aligned
% % and to have the same variance in each dimension as the
% % data.
% sigma = zeros(d,d,M);
% for k = 1:M
%     sigma(:,:,k) = diag(C);
% end

% % Initialize the priors for each of the models
% priors = ones(M,1)/M;

% % Initialize the conditional probabilities for each
% % data point
% % e_z(k,n) = p(x_n|k)
% e_z = ones(M,N)/M;

% for n = 1:300
%     % E Step
%     % ======
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     % Compute <z_ik> as follows:    %
%     %                               %
%     %             p(x_i|k)p(k)      %
%     % <z_ik> = ------------------   %
%     %          sum_j p(x_i|j)p(j)   %
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%     % Compute p(x_i|k)
%     p_x_k = compute_conditional_probability(X,mu,sigma);
%     % Compute p(x_i|k)p(k)
%     joint = p_x_k .* repmat(priors,1,N);
%     % Compute <z_ik>
%     e_z = joint ./ repmat(sum(joint),M,1);

%     % M Step
%     % ======
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     % Update formulae for p(k), mu_k, sigma_k       %
%     %                                               %
%     %        sum_i <z_ik>                           %
%     % p(k) = ------------                           %
%     %             N                                 %
%     %                                               %
%     %        sum_i <z_ik> x_i                       %
%     % mu_k = ----------------                       %
%     %          sum_i <z_ik>                         %
%     %                                               %
%     %           sum_i <z_ik> (x_i-mu_k)(x_i-mu_k)'  %
%     % sigma_k = ----------------------------------  %
%     %                     sum_i <z_ik>              %
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%     % Compute sum_i <z_ik> for all k
%     sum_e_z = sum(e_z')';
%     % Compute the priors
%     priors = sum_e_z/N;
%     % Compute each sigma_k
%     for k = 1:M
%         centered_X = X - repmat(mu(:,k),1,N);
%         sigma(:,:,k) = (centered_X.*repmat(e_z(k,:),d,1))*centered_X'/sum_e_z(k);
%     end
%     % Compute mu_k for all k
%     mu = (X * e_z')./repmat(sum_e_z',d,1);
    
%     figure(3);clf;hold on;
%     for k = 1:M
%         draw_kahal(inv(sigma(:,:,k)),mu(:,k),2);
%     end
% end

% function p_x_k = compute_conditional_probability(X, mu, sigma)
% [d,M] = size(mu);
% [d,N] = size(X);

% p_x_k = zeros(M,N);
% factor = (2*pi)^(-d/2);
% for k = 1:M
%     centered_X = X - repmat(mu(:,k),1,N);
%     p_x_k(k,:) = (factor / det(sigma(:,:,k))^0.5) * exp(-0.5*sum((inv(sigma(:,:,k))*centered_X).*centered_X));
% end
% \end{verbatim}
% \end{appendix}


\end{document}
