\documentstyle[ml99,times]{article}

\def\topfraction{1}
\def\textfraction{0}
\def\dbltopfraction{1}
\def\dbltextfraction{0}
\setcounter{topnumber}{1}
\setcounter{dbltopnumber}{1}

\def\fnum@figure*{\footnotesize {\bf Figure \thefigure}}%%*%%
\def\fnum@figure{\footnotesize {\bf Figure \thefigure}}%%*%%
\def\fnum@table*{\footnotesize {\bf Table \thetable}}%%*%%
\def\fnum@table{\footnotesize {\bf Table \thetable}}%%*%%

\parskip .4pc 
\leftmargini 1.5em
\def\baselinestretch{0.98}

\typeout{############ RUN LATEX 2E ############}

\title{Monte Carlo Hidden Markov Models:\\Learning Non-Parametric
Models of Partially Observable Stochastic Processes}

\author{Sebastian Thrun\ \ \  \ \ \  \ \ \  John C.~Langford\ \ \  \ \ \  \ \ \ Dieter Fox\\[3mm]
Robot Learning Laboratory, School of Computer Science\\
Carnegie Mellon University\\
5000 Forbes Ave., Pittsburgh, PA 15213\\
\{thrun,jcl,dfox\}@cs.cmu.edu}


\begin{document}

\maketitle


\begin{abstract}
\noindent
We present a learning algorithm for non-parametric hidden Markov
models with continuous state and observation spaces.  All necessary
probability densities are approximated using samples, along with
density trees generated from such samples.  A Monte Carlo version of
Baum-Welch (EM) is employed to learn models from data.  Regularization
during learning is achieved using an exponential shrinking
technique. The shrinkage factor, which determines the effective
capacity of the learning algorithm, is annealed down over multiple
iterations of Baum-Welch, and early stopping is applied to select the
right model. Once trained, Monte Carlo HMMs can be run in an any-time
fashion. We prove that under mild assumptions, Monte Carlo Hidden
Markov Models converge to a local maximum in likelihood space, just
like conventional HMMs. In addition, we provide empirical results
obtained in a gesture recognition domain.
\end{abstract}


\section{Introduction}

Hidden Markov models (HMMs) \cite{Rabiner89aABBREV} have been applied
successfully to a large range of temporal signal processing problems,
for tasks such as modeling, filtering, classification and prediction
of partially observable time sequences, stochastic environments.
Hidden Markov models are often the method of choice in areas such as
speech recognition \cite{Rabiner89aABBREV,Waibel90aABBREV}, natural
language processing \cite{Charniak93aABBREV}, robotics
\cite{Meila96aABBREV,Shatkay97aABBREV}, biological sequence analysis
\cite{Henderson96aABBREV,Krogh94aABBREV}, and time series analysis
\cite{Hannaford91aABBREV}.

We present a new hidden Markov model, called {\em Monte-Carlo Hidden
Markov Models (MCHMMs)}.  MCHMMs use samples and non-parametric
density trees to represent probability distributions.  This makes it
possible to learn non-parametric hidden Markov models with continuous
state and observation spaces.  Since continuous state spaces are
sufficiently rich to overfit any data set, our approach uses shrinkage
as a mechanism for regularization. The shrinkage factor, which
determines the effective capacity of the HMM, is annealed down over
multiple iterations of EM, and early stopping is applied for model
selection.

MCHMMs possess the following four characteristics, which
distinguishes them from the traditional HMM approach
\cite{Rabiner89aABBREV}:
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item
{\bf Real-valued spaces.} Both the state space and the observation
space of MCHMMs are continuous. This is important in domains where the
true state and observation space of the environment is continuous.
The importance of continuous-valued spaces has been recognized by
several authors, which have proposed real-valued extensions using
parametric models such as Gaussians and neural networks
\cite{Briegel98aABBREV,Ghahramani98aABBREV,Gordon93aABBREV,Juang85aABBREV,Liporace82aABBREV}.
\item
{\bf Non-parametric.}  Most existing HMM models rely on parametric
densities, defined by a small set of parameters (discrete
distributions included). This is clearly appropriate when one has
reason to believe that the parametric model is correct. In practice,
however, one might lack such strong prior knowledge. MCHMMs learn
non-parametric models of partially observable stochastic systems. Such
non-parametric methods can fit much larger sets of functions, hence
provide additional flexibility.
\item
{\bf Built-in complexity control.} MCHMMs make it possible to control
the complexity of the state space during learning.  In contrast, the
complexity of discrete HMMs is determined by the number of states,
which has to be determined in advance, before training. Multiple
training runs may thus be necessary to determine the right number of
states.
\item
{\bf Resource adaptive.} Finally, after training MCHMMs can be used in
an {\em any-time} fashion. Any-time algorithms can generate an estimate
at {\em any} time; the more time available, however, the better the
answer \cite{Dean88aABBREV,Zilberstein95aABBREV}. The any-time
property is directly inherited from sampling-based methods, as the
number of samples can be controlled on-line.
% \item
% {\bf Computation is placed where needed.} Finally, MCHMM place
% computation ``where needed.'' This property is also inherited from
% sampling-based methods. MCHMM sample in a likelihood-weighted fashion:
% the more likely a state or an observation, the higher the probability
% that it is sampled (and hence requires computation).
\end{enumerate}
This paper presents MCHMMs. It describes a sample-based learning
algorithm that extends Baum-Welch to MCHMMs. We also give theoretical
results that prove (with probability one) convergence of the learning
algorithm to a local maximum in likelihood space. This result, whose
proof is only outlined here (see \cite{Thrun98nABBREV} for all
proofs), extends the well-known convergence result for discrete HMMs
to the real-valued, non-parametric case.  We also present empirical
results obtained in a gesture recognition domain, which seek to
provide the first evidence of the usefulness of MCHMMs in practice.



\section{Generalized Hidden Markov Models}
\label{Sect2}

This section introduces the necessary notation to extend HMMs to
continuous state and observation spaces. It follows the excellent
deposition of \cite{Rabiner89aABBREV}, with continuous density
replacing discrete distributions throughout.  For convergence, we
assume all event spaces and random variables (e.g., state,
observations) are measurable. Unless otherwise specified, all
probability density functions are Lipschitz (implying continuity) and
non-zero over some compact interval.

A {\em generalized HMM} (in short: GHMM) is a partially observable,
time-invariant Markov chain with continuous state and observation
spaces and discrete time.  At each time $t\geq 1$, the HMM's state is
denoted $x_t$.  Initially, at time $t=1$, the state of the HMM is
selected randomly according to the density $\pi$.  State transitions
are governed by a conditional probability density, denoted $\mu(x'\mid
x)$ and called {\em state transition density}.

In HMMs (and thus in GHMMs), the state cannot be observed
directly. Instead, only a probabilistic projection of the state is
observable.  Let $b_{t}$ denote a measurable random variable that
models the observation at time $t$. Observations are generated
according to a probability density conditioned on the state of the HMM
(called the {\em observation density}), denoted $\nu(b\mid x)$ .
A generalized HMM is uniquely defined trough three densities:
$\lambda= \{\pi,\mu,\nu\}$.

Putting computational considerations aside for the moment, knowledge
of $\lambda$ is sufficient to tackle a variety of interesting
practical problems, such as
computing distributions over states and observations at arbitrary
points in time,
generating representative example trajectories in state and
observation space,
determining the
likelihood of example trajectories under an HMM, and
classifying data sequences generated by mixtures of labeled HMMs.
Algorithms for these problems are described in detail in
\cite{Rabiner89aABBREV}; they are easily transferred from the finite
to the continuous case.

In practice, the densities $\lambda$ are often unknown and have to be
estimated from data. The data, denoted $d$, is a sequence of
observations\footnote{For simplicity of the presentation, we only
present the case in which the data consist of a single sequence. The
extension to multiple sequences is straightforward but requires
additional notation.}, denoted
\begin{eqnarray}
d&=&\{O_{1},O_{2}, \ldots, O_{T}\}.
\end{eqnarray}
Here $O_t$ denotes the observation at time $t$. The total number of
observations in $d$ is $T$.

The well-known Baum-Welch algorithm
\cite{Baum66aABBREV,McLachlan97aABBREV,Rabiner89aABBREV} provides a
computationally efficient and elegant approach for learning $\pi$,
$\mu$, and $\nu$.  Baum-Welch begins with an initial model, denoted
$\lambda^{(0)}$. It iterates two steps, an E-step and an M-step (see
also \cite{Dempster77aABBREV}).  In the $n$-th E-step, distributions
for the various state variables $x_t$ are computed under a fixed model
$\lambda^{(n)}$ (with $n\geq 0$). The $n$-th M-step uses these
distributions to derive a new, improved model $\lambda^{(n+1)}$.


In the E-step, distributions are computed for the state variables
$x_t$ conditioned on a fixed model $\lambda$ and the data $d$.  Recall
from \cite{Rabiner89aABBREV} that in the discrete case,
\begin{eqnarray}
\alpha_t^{(n)}(x) &=& Pr(x_t=x\mid O_1,\ldots,O_t,\lambda^{(n)})\\
\beta_t^{(n)}(x) &=& Pr( O_{t+1},\ldots,O_{T}\mid x_t=x, \lambda^{(n)})\\
\gamma_t^{(n)}(x) &=& Pr(x_t=x \mid d, \lambda^{(n)})\\
\xi_t^{(n)}(x,x') &=& Pr(x_t=x , x_{t+1}=x' \mid d, \lambda^{(n)})
\end{eqnarray}
The continuous case is analogous; however, here $\alpha$ and $\beta$
are densities (and thus may be larger than 1).  Following
\cite{Rabiner89aABBREV}, these densities are computed incrementally;
$\alpha$ is computed forward in time, and $\beta$ backwards in time
(for which reason this algorithm is often referred to as the {\em
forward-backward algorithm}). Initially,
\begin{eqnarray}
\alpha_0^{(n)}(x)\;\;=\;\;\pi^{(n)}(x)\;\; &\mbox{and}&\;\;  \beta_{T}^{(n)}(x)\;\;=\;\; 1
\label{alphaInit}
\label{betaInit}
\end{eqnarray}
and for all other $\alpha_t$ and $\beta_t$:
\begin{eqnarray}
\hskip -12mm
&&\alpha_t^{(n)}(x)= \int \alpha_{t-1}^{(n)}(x') \;\mu^{(n)}(x\mid x')\; \nu^{(n)}(O_t \mid x) \;dx'
\label{eqAlpha}\\
\hskip -12mm
&&\beta_t^{(n)}(x) = \int \beta_{t+1}^{(n)}(x') \;\mu^{(n)}(x'\mid x)\;\nu^{(n)}(O_{t+1} \mid x')\;dx'
\label{eqBeta}
\end{eqnarray}
Bayes rule governs the conditional density over the state space at
time $t$:
\begin{eqnarray}
\gamma_t^{(n)}(x)&=&\frac{\alpha_t^{(n)}(x)\;\beta_t^{(n)}(x)}{\displaystyle \int \alpha_t^{(n)}(x')\;\beta_t^{(n)}(x')\;dx'}
\label{eqGamma}
\end{eqnarray}
Similarly, the state transition densities $\xi^{(n)}$ are computed as
\begin{eqnarray}
\xi_t^{(n)}(x,x') &=&\\
&&\hskip -25mm \frac{\alpha_{t}^{(n)}(x)\;\mu^{(n)}(x'\mid x)\;\nu^{(n)}(O_{t+1}\mid x')\;\beta_{t+1}^{(n)}(x')}{\displaystyle \int\int \alpha_{t}^{(n)}(\bar{x})\;\mu^{(n)}(\bar{x}'\mid \bar{x})\;\nu^{(n)}(O_{t+1}\mid \bar{x}')\;\beta_{t+1}^{(n)}(\bar{x}')\;d\bar{x}\;d\bar{x}' }
\nonumber
\end{eqnarray}
This computation is completely analogous to the finite case, replacing
conditional probabilities by conditional densities.

The M-step uses $\gamma_t^{(n)}(x)$ and $\xi_t^{(n)}(x,x')$ to compute
a new model $\lambda^{(n+1)}$, using the maximum likelihood estimator:
\begin{eqnarray}
\pi^{(n+1)}(x) &=& \gamma_0^{(n)}(x)
\label{eqPi}\\
\mu^{(n+1)}(x'\mid x) &=& \frac{ \sum_{t=1}^{T-1} \xi_t^{(n)}(x,x')}{ \sum_{t=1}^{T-1} \gamma_t^{(n)}(x)}
\label{eqMu}\\
\nu^{(n+1)}(b\mid x) &=&  \frac{ \sum_{t=1}^{T} I_{O_t=b} \; \gamma_t^{(n)}(x)}{ \sum_{t=1}^{T} \gamma_t^{(n)}(x)}
\label{eqNu}
\end{eqnarray}
Here $I_{cond}$ denotes an indicator variable that is 1 if the
condition {\em cond} is true, and 0 otherwise.  A straightforward
result is the convergence of GHMMs.

{\bf Theorem~1. (GHMM Convergence Theorem)} {\sl Under the assumptions
above, the GHMM learning algorithm does not decrease the likelihood
of the observation sequence with probability 1. With probability 1, it
does not improve the likelihood if, and only if a local maximum or a
saddle point has been reached.}

The full proof can be found in \cite{Thrun98nABBREV}. Our proof is a
direct extension of a result by Juang \cite{Juang85aABBREV}, who has
shown convergence to local maxima for HMMs with a finite number of
states and continuous observations, where the observation densities
are a mixture of log concave or ellipsoidal symmetrical densities. The
convergence of GHMMs follows from the fact that any continuous and
Lipschitz density can be approximated by mixtures of Gaussians, hence
by a sequence of HMMs that meet Juang's assumptions. 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Sampling and Density Trees}

\subsection{Sampling}
\label{SectSamp1}

MCHMMs use samples to approximate all densities (models, posteriors,
conditional and joint densities).  Sample sets are (finite) sets of
values, annotated by numerical probabilities
\cite{Liu98aABBREV,Neal93aABBREV}.  More specifically, let $f$ be a
probability density function, and let $N$ denote a positive number
(the cardinality of a sample set).  A {\em sample set} is a
set
\begin{eqnarray}
X&=&\{\langle x_{1},p_{x_1}\rangle\, \ldots,\langle x_{N},p_{x_N}\rangle \}
\end{eqnarray}
where for all $n$ with $1\leq n \leq N$: $x_{n}\in {\rm dom}(f)$,
$p_{x_n}\in [0,1]$, and $\sum_{n=1}^N p_{x_n}=1$.  Sample sets can be
thought of as discrete probability distributions over the event space
$\{x_{1},\ldots,x_{N}\}$, defined through the $p$-values
\cite{Neal93aABBREV,Pitt99aABBREV} For reasons that will become clear
below, $p$-values are often referred to as {\em importance factors}
\cite{Rubin88aABBREV}.

A {\em sampling method} is a method that approximates $f$ through
samples.  We will call a sampling method {\em asymptotically
consistent} if for $N \rightarrow \infty$, $X$ converges to $f$ with
probability 1 when integrated over the system of half-open Borel sets:
\begin{eqnarray}
\lim_{N \rightarrow \infty}\sum_{\langle x,p_x\rangle\in X} I_{x_0\leq x < x_1} \; p_x
&= &
\int_{x_0}^{x_1} f(x)\;dx
\label{eqConv}
\end{eqnarray}
Sampling can be applied to approximate any probability density $f$
\cite{Tanner93aABBREV}.




%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.1mm}
\begin{figure*}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,67)(0,0)
\put(1,17){\special{psfile=pictures/hmm-21.ps vscale=25.3 hscale=25.3}}
\put(1,-17){\special{psfile=pictures/hmm-20.ps vscale=25.3 hscale=25.3}}
\put(45,-40){\special{psfile=pictures2/hmm.1.pdf.rot.ps vscale=41.8 hscale=41.8}}
\put(3,63){\footnotesize (a)}
\put(3,29){\footnotesize (b)}
\put(55,60){\footnotesize (c)}
\end{picture}
\caption{(a) Data set (b) partitioned by a density tree.  (c) Density
tree with shrinkage for complexity control. Each branch cuts the area
under a node into two equally-sized, rectangular halves, shown in
grey. The density of each node is an exponentially weighted mixture of
densities of the nodes along the path to the leaf. The $\rho$-terms
describe the mixture coefficients for this example.}
%
\label{FigSamplesAndTree}%key%
%\rule{\textwidth}{.2mm}
\ \\[-10mm]
\end{figure*}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%

We will distinguish two methods for sampling, both of which are
asymptotically consistent:
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item
{\bf Likelihood-weighted sampling.}  When $f$ is available explicitly,
we will directly sample from $f$ (e.g., using rejection sampling
\cite{Tanner93aABBREV}).  In likelihood-weighted sampling, all
$p$-values are $1/N$. The concentration of the resulting sample set
will be proportional to the density $f$, i.e., the larger the measure
of $f$ over a region $A$, the more samples will be in $A$ (in
expectation).
\item
{\bf Importance sampling.}  In some cases, however,
likelihood-weighted sampling cannot be applied. This is the case, for
example, when approximating $\alpha_t^{(n)}(x)$, which is a recursive
product of convolved densities and other densities conditioned on
observations (c.f., (\ref{eqAlpha})). Instead of sampling directly
from $f$, importance sampling \cite{Rubin88aABBREV} samples from a
different distribution $g$ that has the property $f(x)>0
\Longrightarrow g(x)>0$. The difference between $f$ and $g$ is
accounted by the $p$-values, which are defined as follows:
\begin{eqnarray}
\bar{p}_{x_n}&=& f(x_n)/g(x_n)
\\
p_{x_n}&=& \frac{\bar{p}_{x_n}}{\sum_{n=1}^N \bar{p}_{x_n}}
\end{eqnarray}
Importance sampling \cite{Rubin88aABBREV} has been studied extensively
in the statistics literature (see for example \cite{Tanner93aABBREV}).
The reader may notice that likelihood-weighted sampling is a special
case of importance sampling (with $f\equiv g$).
\end{enumerate}
Figure~\ref{FigSamplesAndTree}a shows a sample set drawn by
likelihood-weighted sampling from a distribution that resembles the
shape of a sine wave. All probabilities $p_x$ of the sample set
shown there are the same, and the samples are concentrated in a small
region of the $\Re^2$.

Both sampling methods can equally be applied to {\em sample from a
sample set} (often referred to as {\em resampling}
\cite{Rubin88aABBREV}).  Let $X$ be a sample set. Likelihood-weighted
sampling draws samples from $X$ according to the (discrete)
probability distribution induced by their $p_x$-values. Each sample is
then assigned the same probability ($1/N$).  Importance sampling draws
samples $x$ from $X$ according to some other density $g$, and assigns
$\frac{p_x}{g(x)}$ to the new sample's $p$-value.  Sampling from
sample sets plays an important role in the Monte Carlo HMM algorithm
described below.

A straightforward but important result is the following: 

{\bf Theorem~2.} {\sl Likelihood-weighted sampling converges (in the
sense of (\ref{eqConv}) with probability one to $f$ as
$N\longrightarrow \infty$, at a rate of $1/\sqrt{N}$. Importance
sampling also converges with probability one to $f$ at the same
asymptotic rate if $\frac{f(x)}{g(x)} < \infty$.}

The proof, which is given in detail for likelihood-weighted sampling
in \cite{Thrun98nABBREV}, follows directly from a result in
\cite{Tanner93aABBREV}, pg.~33.  For the importance sampler, the
variance of the error depends on the ``mismatch'' between $f$ and $g$.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





\subsection{Particle Filters}
\label{SectPartFil}

Armed with the definition of sample sets and procedures for generating
them, we can now devise a sample-based algorithm for tracking the
state of a dynamical system: {\em particle filters}
\cite{Pitt99aABBREV}.  Variants of this algorithm, which have become
popular in recent years, can be found in the literature under names
such as {\em survival of the fittest} \cite{Kanazawa95aABBREV},
{\em condensation algorithm} \cite{Isard98aABBREV}, or, in the context of
mobile robotics, {\em Monte Carlo localization}
\cite{Dellaert99aABBREV,Fox99aABBREV}.

Particle filters assume that the model $\lambda$ is known. They
approximate the distributions of $\alpha_t$ using sampling/importance
resampling:
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item {\bf Initialization:} Generate $N$ samples of $\pi$ (e.g., using
likelihood-weighted sampling).
\item {\bf For each $t=1,\ldots $ do:}
\begin{enumerate}\itemsep 0mm\parskip 0mm \parindent 0mm
\item
  Generate $N$ samples $\langle x',p_{x'}\rangle$ from the sample set
  representing $\alpha_{t-1}$ using likelihood-weighted sampling.
\item
  For each sample $\langle x',p_{x'}\rangle$, generate a random $x$
  from $\mu(x\mid x')$ using likelihood-weighted sampling.
\item
  Set $p_{x}$ to a value proportional to $\nu(O_t\mid x)$, where $O_t$
  is the $t$-th observation in the data set. 
\end{enumerate}
\end{enumerate}
Particle filters have successfully been applied to state estimation
problems in computer vision \cite{Isard98aABBREV} and robotics
\cite{Dellaert99aABBREV,Fox99aABBREV}.  In MCHMMs, this algorithm is
used to approximate the distributions of $\alpha_t$ and $\beta_t$.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Density Trees}
\label{SectDensTrees}

While sample sets are sufficient to approximate con\-tinuous-valued
distributions, they differ from those in that they are {\em
discrete}. This is problematic if one wants to combine
densities represented through different sample sets. For example,
according to Equation (\ref{eqGamma}), $\gamma_t$ is the normalized
product of two densities, $\alpha_t$ and $\beta_t$, each of which
represented by its own sample set. Unfortunately, with probability
one, none of the samples in $\alpha_t$ and $\beta_t$ are identical,
and thus it is {\em not} straightforward how to obtain an
approximation of their product.

MCHMMs solve this problem by transforming sample sets into
(non-parametric) density trees
\cite{Koller98aABBREV,Moore97aABBREV,Omohundro91aABBREV}, which
effectively ``generalize'' discrete distributions (samples) to
continuous distributions.  Each node in a density tree is annotated
with a hyper-rectangular subspace of ${\rm dom}(f)$, denoted $V$ (or
$V_i$ for the $i$-th node).  Initially, all samples are assigned to
the root node, which covers the entire domain of $f$.  The tree is
grown by splitting each node $i$ whenever the following two conditions
are fulfilled:
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item 
At least $\sqrt{N}$ samples $\langle x,p_x\rangle\in X$ fall into 
$V_i$.
\item Its depth, i.e., its distance from the root node, does not
exceed $\lfloor \frac{1}{4} \log_2 N \rfloor$.%\\[-7mm]
\end{enumerate}
If a node is split, its hyper-rectangle $v$ is divided into two
equally sized hyper-rectangles along its longest dimension, as
illustrated in
Figure~\ref{FigSamplesAndTree}c. Figure~\ref{FigSamplesAndTree}b shows
a density tree generated for the sample set shown in
Figure~\ref{FigSamplesAndTree}a.


Let $f$ be a density function, $X$ a sample drawn from this density,
and for all $x\in {\rm dom}(f)$ let $i(x)$ denote the leaf whose
region $V_i$ contains $x$.  Furthermore, let $\hat{\sigma}_i$ denote
the relative frequency of samples in the $i$-th node weighted by their
respective $p_x$ values:
\begin{eqnarray}
\hat{\sigma}_i &=& \sum_{\langle x,p_{x}\rangle\in X}I_{x\in V_{i}}\; p_{x}
\label{eqSigma}
\end{eqnarray}
Hence the density function of a tree, denoted $\hat{f}$, is defined as
follows:
\begin{eqnarray}
\hat{f}(x) &=& \frac{\hat{\sigma}_{i(x)}}{| V_{i(x)} |}\;\;\;\;\;\;\;\;\;\;(\forall x\in {\rm dom}(f))
\label{eqTreeDens}
\end{eqnarray}
The numerator of (\ref{eqTreeDens}) describes the weighted relative
frequency that a sample falls into the interval of the node
$i(x)$. The denominator is the size of the interval, which is
necessary to transform a relative frequency into a density over the
interval $V_{i(x)}$.  The density function $\hat{f}(x)$ can be equally
defined for {\em internal nodes} (non-leaf nodes). This will be
important below, where estimates at different levels of the tree are
mixed for regularization.


In \cite{Thrun98nABBREV} we prove the following important convergence
result:

{\bf Theorem~3.} {\sl With probability one, density trees converge to
$f$ as $N\mbox{$\rightarrow$} \infty$.}

The proof exploits the fact that our splitting rule ensures that (1)
the depth of the tree grows without bounds and (2) the number of
samples in each leaf also grows without bound as $N\mbox{$\rightarrow$}
\infty$. The convergence then follows from the fact that $f$ is
Lipschitz.\footnote{The termination conditions for growing trees where
chosen to facilitate the derivation of Theorem~3. In practice, these
conditions can be overly restrictive, as they often require large
sample sets to grow reasonable-sized trees. Our actual implementation
sidesteps these conditions, and trees are grown all the way to the
end. While the convergence results reported here are not applicable
any longer, our implementation yielded much better performance
specifically when small sample sets were used (e.g., 100 samples).}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%#@a-------------------------------------------------------------------------%
\begin{table*}[t]
\begin{footnotesize}
\begin{description}
\item[Model Initialization:]

Initialize $\lambda= \{\pi,\mu,\nu\}$ by three randomly drawn sets of
samples of the appropriate dimension. Generate density trees from
these samples. Set $\rho=1$, and chose an initial sample
set size $N > 0$.

\item[E-step:]
\ \\[-4mm]
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item
 Computation of $\alpha_0$ (c.f., (\ref{alphaInit})). Generate $N$
 samples from the tree representing $\pi$ using likelihood-weighted
 sampling.
\item
Computation of $\alpha_t$ (c.f., (\ref{eqAlpha})). For each $t$ with
$1<t \leq T$ do:
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item
  Generate $N$ samples $\langle x',p_{x'}\rangle$ from the sample set
  representing $\alpha_{t-1}$ using likelihood-weighted sampling.
\item
  For each sample $\langle x',p_{x'}\rangle$, generate the conditional
  density $\mu(x\mid x')$ using the tree-version of $\mu$.  Sample a
  single $x$ from this tree, using likelihood-weighted sampling.
\item
  Set $p_{x}$ to a value proportional to $\nu(O_t\mid x)$, where $O_t$
  is the $t$-th observation in the data set.  This density value is
  obtained using the tree representing $\nu$.
\item
  Generate a tree from the new sample set.

\end{enumerate}
\item
 Computation of $\beta_T$ (c.f., (\ref{betaInit})). Generate $N$
 uniformly distributed samples.
\item
  Computation of $\beta_t$ (c.f., (\ref{eqBeta})). For each $t$ with
  $T>t \geq 1$ do:
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item
  Generate $N$ samples $\langle x',p_{x'}\rangle$ from the sample set
  representing $\beta_{t+1}$ using likelihood-weighted sampling.
\item
  For each sample $\langle x',p_{x'}\rangle$, generate the conditional
  density $\mu(x'\mid x)$ using the tree-version of $\mu$.  Sample a
  single $x$ from this tree, using likelihood-weighted sampling.
\item
  Set $p_{x}$ to a value proportional $\nu(O_{t+1}\mid x')$, where
  $O_{t+1}$ is the $t+1$-th observation in the data set.  This density
  value is obtained using the tree representing $\nu$.
\item
  Generate a tree from the new sample set.
\end{enumerate}

\item
  Computation of $\gamma_t$ (c.f., (\ref{eqGamma})). For each $t$ with
  $1\leq t \leq T$ do:
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item
  Generate $N/2$ sample from $\alpha_t$ by likelihood weighted
  sampling, and assign to each sample $x$ an importance factor $p_x$
  proportional to $\beta_t(x)$, using the tree approximation of
  $\beta_t$.
\item 
  Generate $N/2$ sample from $\beta_t$ by likelihood weighted sampling
  and assign to each sample $x$ an importance factor $p_x$
  proportional to $\alpha_t(x)$, using the tree approximation of
  $\alpha_t$.
\end{enumerate}
\end{enumerate}
\ \\[-12mm]
\item[M-step:]
\ \\[-4mm]
\begin{enumerate}%\itemsep 0mm\parskip 0mm \parindent 0mm
\item
    Estimation of the new state transition density $\mu$ (c.f.,
    (\ref{eqMu})): Pick $N$ random times $t\in\{1,\ldots,T-1\}$ and
    generate samples $\langle x,p_x \rangle $ and $\langle x',p_{x'}
    \rangle $ from $\gamma_t$, and $\gamma_{t+1}$, respectively, by
    likelihood-weighted sampling. Add $\langle (x,x'),N^{-1} \rangle $
    into the sample set representing $\mu$.  Generate a tree from the
    sample set. When using $\mu(x'\mid x)$, condition the tree on $x$.
\item    
    Estimation of the new observation density $\nu$ (c.f.,
    (\ref{eqNu})): Pick $N$ random $t\in\{1,\ldots,T\}$ and
    generate a sample $\langle x,p_x \rangle $ from $\gamma_t$ by
    likelihood-weighted sampling. Add $\langle (x,O_t),N^{-1} \rangle $
    into the sample set representing $\nu$. Generate a tree from the
    sample set. When using $\nu(b\mid x)$, condition the tree on $x$.
\item
    Estimation of the new initial state distribution $\pi$ (c.f.,
    (\ref{eqPi})): Copy the sample set $\gamma_0$ into $\pi$. Generate
    a tree from the sample set.
\end{enumerate}
\item[Annealing:] Set $\rho \leftarrow \rho \bar{\rho}$. Stop when the
likelihood of an independent cross-validation set is at its maximum.
\item[Sample set size:] $N \longleftarrow \eta N$, for some $\eta > 1$.
\end{description}
\end{footnotesize}
\ \\[-13mm]
\caption{The MCHMM algorithm.}
\label{TabMCHMM}
%\rule{\textwidth}{.2mm} 
\ \\[-10mm]
\end{table*}
%#@b-------------------------------------------------------------------------%


\subsection{Regularization Through Shrinkage and Annealing}

The continuous nature of the state space in GHMMs (and hence in
MCHMMs), makes it possible to overfit {\em any} data set, no matter
how large. This is because GHMMs are rich enough to assign a different
state to each observation in the training data (of which there are
only finitely many), making it essentially impossible to generalize
beyond sequences other than the ones presented during training.  A
similar problem arises in conventional HMMs, if they are given more
states than there are training examples. In GHMMs the problem is
inherent, due to the topology of continuous spaces.  Thus, some kind
of regularization is needed to prevent overfitting from happening.

Our approach to regularization is based on {\em shrinkage}
\cite{Stein55aABBREV}.  Shrinkage is a well-known statistical
technique for ``lumping together'' different estimates from different
data sources. In a remarkable result by Stein \cite{Stein55aABBREV},
shrinking estimators were proven to yield uniformly better solutions
over unbiased maximum-likelihood estimators in multivariate Gaussian
estimation problems. {\em Shrinkage trees} were introduced by
McCallum \cite{McCallum98aABBREV}.  Instead of using the density
estimates at the leafs of a tree, shrinkage trees mix those densities
with densities obtained further up in the tree. These internal
densities are less specific to the region covered by a leaf node;
however, their estimates are usually obtained from more data, making
them less susceptible to variance in the data.

Figure~\ref{FigSamplesAndTree}c shows an example using shrinkage with
an exponential factor, parameterized by $\rho$ (with $0\leq\rho\leq
1$). Here each node in the tree, with the exception of the root node,
weighs its own density estimate by $(1-\rho)$, and mixes it with the
density estimate from its parent using the weighting factor $\rho$. As
a result, every node along the path contributes to the density
estimate at the leaf; however, its influence decays exponentially with
the distance from the leaf node.  Obviously, the value of $\rho$
determines the amount of shrinkage. If $\rho=1$, only the root node is
consulted, hence, the probability density induced by the tree is
uniform over its domain. If $\rho=0$, on the other hand, there is no
shrinkage and only the estimates in the leaf nodes determine the shape
of the density function. For intermediate values of $\rho$, estimates
along the entire path are combined.

Since the optimal value of $\rho$ is problem-specific---de facto it
depends on the nature of the (unobservable) state space---our approach
uses annealing and cross validation to determine the best value for
$\rho$. More specifically,
\begin{eqnarray}
\rho^{(n)} &=& \bar{\rho}^{\;n-1}
\end{eqnarray}
where $\bar{\rho} < 1$ is a constant (e.g., 0.9) and $n$ denotes the
iteration counter of the Baum-Welch algorithm (starting at 1). Thus,
$\rho$ starts with $\rho^{(0)}=1$, for which nothing can be learned,
since every density is uniform. The parameter $\rho$ is then annealed
towards zero in an exponential fashion.  Cross validation (early
stopping) is applied to determine when to stop training.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.0mm}
\begin{figure*}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,36)(25,33)
\put(-6,-20){\special{psfile=pictures/gest-20.ps vscale=50 hscale=50}}
\put(36,-20){\special{psfile=pictures/gest-21.ps vscale=50 hscale=50}}
\put(78,-20){\special{psfile=pictures/gest-22.ps vscale=50 hscale=50}}
\put(120,-20){\special{psfile=pictures/gest-23.ps vscale=50 hscale=50}}
\end{picture}
\caption{Examples of gestures. The first two gestures were drawn
counterclockwise; whereas the other two were drawn clockwise. The
MCHMM learning task is to differentiate them.}
%
\label{FigData}%key%
%\rule{\textwidth}{.2mm}
\ \\[-10mm]
\end{figure*}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%

%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.11111mm}
\begin{figure*}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,48)(-8,2)
\put(-18,-41){\special{psfile=pictures/hmm-1c-mod.ps vscale=52 hscale=50}}
\put(60,-41){\special{psfile=pictures/hmm-1-mod.ps vscale=52 hscale=50}}
\put(1,40){\footnotesize (a)}
\put(79,40){\footnotesize (b)}
\end{picture}
\caption{Log likelihood  as a function of the iteration of EM, for the
synthetic data set.  {\bf (a) Training}: Top curve: 1000 samples for
all densities. Middle curve: 1000 samples for $\mu$ and $\nu$, 100
samples for $\alpha$, $\beta$, $\gamma$, and $\pi$.  Bottom curve: 100
samples for $\mu$ and $\nu$, 10 samples for $\alpha$, $\beta$,
$\gamma$, and $\pi$.  {\bf (b) Testing}: The top curve shows the log
likelihood of ``positive'' training data (same orientation as used for
training), whereas the bottom curve shows the log likelihood for the
``negative'' class (opposite orientation as used for training).  Each
result is averaged over 10 independent experiments; 95\% confidence
bars are also shown.}
%
\label{FigRes0}%key%
%\rule{\textwidth}{.2mm}
%\ \\[-13mm]
%\end{figure*}
%\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%
\bigskip
%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.11111mm}
%\begin{figure*}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,48)(-8,0)
\put(-18,-41){\special{psfile=pictures/hmm-3-mod.ps vscale=50 hscale=50}}
\put(60,-41){\special{psfile=pictures/hmm-4-mod.ps vscale=50 hscale=50}}
\put(1,40){\footnotesize (a)}
\put(79,40){\footnotesize (b)}
\end{picture}
\caption{(a) Log likelihood for the gesture data base obtained for
independent testing data. The top curve shows the log likelihood of
gestures drawn in same same way as the training gesture, whereas the
bottom curve shows the log likelihood of gestures drawn in opposite
direction. (b) Trade-off between sample set size and log likelihood, for
the gesture data set. The more samples are available, the higher the
likelihood, and hence the more accurate the results. However, even
extremely small sample sets yield reasonable discrimination.}
%
\label{FigRes2}%key%
%\rule{\textwidth}{.2mm}
\ \\[-10mm]
\end{figure*}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%



\section{Monte Carlo HMMs}
\label{Sect5}
We are now ready to present the main algorithm of this paper, along
with the main theoretical result: The Monte Carlo algorithm for GHMMs,
called {\em Monte Carlo hidden Markov models} (in short {\em
MCHMM}). An MCHMM is a computational instantiation of a GHMM that
represents all densities through samples and trees.  

The learning algorithm for MCHMMs is depicted in Table~\ref{TabMCHMM}.
MCHMMs use both sample-based and tree representations during
learning.  In a nutshell,
\begin{itemize}
\item
{\bf E-step.} Particle filters are applied for forward and backward
projection (c.f., Section \ref{SectPartFil} and Equations
(\ref{eqAlpha}) and (\ref{eqBeta})), using the tree approximations of
$\pi$, $\mu$ and $\nu$.

Forward and backward densities are multiplied (c.f., Equation
(\ref{eqGamma})) by multiplying the importance factors of one of the
sample sets (e.g., $\alpha_t$) with their density values under the
tree induced from the other sample set (e.g., $\beta_t$). This
``trick'' overcomes the problem discussed in the beginning of Section
\ref{SectDensTrees}, namely that with probability 1 no two samples in
the sets representing $\alpha_t$ and $\beta_t$ are the same. Its
asymptotic correctness follows from the convergence results for
sampling and density trees, and the convergence of importance sampling
\cite{Rubin88aABBREV,Tanner93aABBREV}.
\item
{\bf M-step.} The new model $\lambda=\{\pi,\mu,\nu\}$ is obtained by
generating samples from the initial state distribution (in the case of
$\pi$), distributions of pairs of consecutive states (in the case of
$\mu$), and distributions of states paired with the corresponding
observations (in the case of $\nu$). Trees are then induced from these
sample sets. 

The trees of $\mu$ and $\nu$ represent joint distributions, not
conditional distribution. Hence, when using these trees they are
conditioned on a state $x$ (c.f., Equations (\ref{eqMu}) and
(\ref{eqNu})).
\end{itemize}
For regularization, MCHMMs use annealing and cross-validation to
determine the best shrinkage factor.  To ensure convergence, the
number of samples $N$ is increased over time.

The convergence of MCHMMs is guaranteed by the following important
theorem:

{\bf Theorem~4 (MCHMM Convergence Theorem).} {\sl Under the
assumptions stated above, any strictly monotonically increasing
schedule of $ N $'s will cause an MCHMM to converge to a local maximum
of the corresponding GHMM with probability 1. }

We omit the lengthy proof, which can be found in \cite{Thrun98nABBREV}.

Once an MCHMM has been trained, it suffices to memorize only the
tree-based version of the model $\lambda= \{\pi,\mu,\nu\}$; all sample
sets can be discarded. Problems such as state tracking, prediction,
analysis and discrimination can then typically be solved using
particle filters (or variants/extensions thereof). By adapting the
sample set size on-line, these methods can be implemented as {\em
any-time} algorithms \cite{Dean88aABBREV,Zilberstein95aABBREV} that
adapt their computational complexity to the available resources. In
addition, sampling in a likelihood-weighted manner places computation
where needed: at regions with high likelihood.  These properties make
MCHMMs attractive for real-time tracking and control
\cite{Fox99aABBREV}.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%








\section{Experimental Results}
\label{Sect7}
This section reports experimental results, obtained for MCHMMs applied
to two problems: an artificial one which was chosen for demonstration
purposes, and a gesture recognition problem. The results are primarily
intended to illustrate MCHMMs in the finite sample case. Comparisons
with existing HMM methods is beyond the scope of this paper.

The artificial data set consists of multiple oscillatory
sequences. Observations are 10-dimensional and governed by
\begin{eqnarray}
O_t =\left\{\begin{array}{ll}
\langle 0.25+\varepsilon_{t,1}, 0.25+\varepsilon_{t,2},0.25+\varepsilon_{t,3},\\
\ \ \  \ldots,0.25+\varepsilon_{t,10}\rangle &
\mbox{if $t$ odd} \\
& \\
\langle 0.75+\varepsilon_{t,1},0.75+\varepsilon_{t,2},0.75+\varepsilon_{t,3}, \\
\ \ \ \ldots,0.75+\varepsilon_{t,10}\rangle &
\mbox{if $t$ even} 
\end{array}\right.
\nonumber
\end{eqnarray}      
where $\varepsilon_{t,i}$ (with $1\leq t \leq 20$ and $1\leq i \leq
10$) are independent and identically distributed noise variables with
zero-centered triangular distribution over $[-0.15;0.15]$. To test the
discriminatory accuracy of the learned model, we also generated a
second, similar data set (using other independent noise variables
$\varepsilon$):
\begin{eqnarray}
O_t =\left\{\begin{array}{ll}
\langle 0.25+\varepsilon_{t,1}, 0.75+\varepsilon_{t,2},0.25+\varepsilon_{t,3},\\
\ \ \ \ldots, 0.75+\varepsilon_{t,10}\rangle &
\mbox{if $t$ odd} \\
& \\
\langle 0.75+\varepsilon_{t,1},0.25+\varepsilon_{t,25},0.75+\varepsilon_{t,3},\\
\ \ \ \ldots, 0.25+\varepsilon_{t,10}\rangle &
\mbox{if $t$ even} 
\end{array}\right.
\nonumber
\end{eqnarray}      
The second data set consisted of a collection of hand-drawn
gestures. Figure~\ref{FigData} shows examples. Once drawn, all
gestures in our data base look quite similar.  However, some of the
gestures were drawn clockwise, whereas others were drawn
counterclockwise. Here we are interested in discriminating clockwise
from counterclockwise gestures. Notice that this problem is more
difficult then the artificial one, as the observations alone (stripped
of their temporal order) are insufficient for discrimination; instead,
the MCHMM has to learn a meaningful model of the internal state.



Figures~\ref{FigRes0} shows results obtained for the first dataset.
Shown in both figures are curves that characterize the log likelihood
(notice that the log likelihood is unbounded for continuous
distributions). Figure~\ref{FigRes0}a illustrates that the log likelihood of
the {\em training set} increases monotonically over multiple
iterations of EM (Baum-Welch), for different experimental conditions
described in the figure caption.  Figure~\ref{FigRes0}b shows testing
results.  The upper curve depicts the log likelihood for a set of random
sequences that are generated from the same model as the training
sequences.  The bottom curve depicts the log likelihood for independently
generated sequences using the other model, for which the MCHMM was not
trained.  In both cases, only a single data set (of length 20) was
used for training, and testing was performed using 20 data sets
generated independently.  The initial shrinkage factor was $\rho=1$,
which was annealed down at the rate $\bar{\rho}=0.9$.
Figure~\ref{FigRes0}b shows the effect of overfitting: the testing
likelihood increases first, but then levels off (approx.~iteration 6)
and decreases. Using early stopping, we obtained 100\% generalization
accuracy in all experiments, even if as few as $N=10$ samples were
used.


Figure~\ref{FigRes2}a shows the corresponding testing result for the
more difficult gesture recognition problem. These results were
obtained using a single gesture for training only, and using 50
gestures for testing.  The top curve depicts the log likelihood for
gestures drawn in the same direction as the training gesture, whereas
the bottom curve shows the log likelihood for gestures drawn in the opposite
direction.  As easily noticed in Figure~\ref{FigRes2}a, the difference
between classes is smaller than in Figure~\ref{FigRes0}b---which comes
at no surprise---yet the likelihood of the ``correct'' class is still a
factor of 4 to 5 larger than that of the opposite
class. Figure~\ref{FigRes2}a also shows the effect of annealing.  The
best shrinkage value is obtained after 12 iterations, where
$\rho=0.28$. As in the previous case, cross-validation using a single
gesture performs well. On average, the classification accuracy when
using cross-validation for early stopping is 86.0\%. This rate is
remarkably high, given that only a single gesture per class was used
for training.


Figure~\ref{FigRes2}b illustrates the trade-off between computation
and accuracy empirically for the gesture data base.  Shown there is
the log likelihood of the testing data as a function of the sample set
size $N$. Notice that the horizontal axis is logarithmic.  All of
these points are generated using a model $\lambda$ obtained using
early stopping.  The sample set size was generated {\em after}
training, to investigate the effect of computational limitations on
the on-line performance of the MCHMM. As in Figure~\ref{FigRes2}a, the
top curve in Figure~\ref{FigRes2}b depicts the log likelihood of
gestures drawn in the same direction as the training data, whereas the
bottom curve shows the log likelihood of gestures drawn in opposite
direction.  The result in Figure~\ref{FigRes2}b illustrates that the
likelihood (and hence the accuracy) of both data sets increases with
the sample set size.  This is not surprising, as the accuracy of the
approximations increases with the sample set size. In addition,
Figure~\ref{FigRes2}b suggests that the distance between clockwise and
counterclockwise gestures is approximately constant in log likelihood
space (hence it grows in likelihood space). This illustrates that better
results are achieved with larger sample sets. In our experiments,
however, even small sample sets yielded good discrimination (after
training with larger sample sets).  For example, we observed on
average 16\% classification error with $N=10$ samples.

For the gesture data base, one iteration of EM (using the settings
described in the caption of Figure \ref{FigRes0}a) took the following
time:
\begin{center}
\begin{tabular}{c|c|c}
sample set size & sample set size & training time\\[-0.5mm]
for  $\mu, \nu$ &
 for $\alpha, \beta, \gamma, \pi$ &
(in seconds)\\
\hline
$\;\,100$ & $\;\,100$ & 0.852 \\
$1000$ & $\;\,100$ & 1.53 \\
$\;\,100$ & $1000$ & 10.36 \\
$1000$ & $1000$ & 18.38 \\
\end{tabular}
\end{center}
After training, processing a single gesture (which involves running
the particle filter on two different MCHMMs) took the following time:
\begin{center}
\begin{tabular}{c|c}
sample set size for $\alpha$ & running time (in seconds)\\
\hline
$\;\;\,\,10$ & $0.00242$ \\
$\;\,100$ & $0.0242\;\,$ \\
$1000$ & $0.260\;\,\;\,$
\end{tabular}
\end{center}
These numbers were obtained on a low-end PC (equipped with a 266Mhz
Pentium II). The latter numbers illustrate that after learning, the
processing of a gesture is several orders of magnitude faster than the
hand motion involved in gesture generation.


% 10 samples, 0.00134 secs
% 100 samples, 0.023 secs
% 1000 samples, 1.40 secs


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



\section{Conclusion}
\label{Sect9}
We have presented Monte Carlo Hidden Markov Models (MCHMMs), a new
algorithm for hidden Markov model learning. MCHMMs extend HMMs to
real-valued state and observation spaces and non-parametric
models. They represent all densities by samples, which are transformed
into probability density functions using density trees.  Because the
continuous state spaces are rich enough to fit (and overfit!)
arbitrary data sets, our approach uses shrinkage to reduce its
complexity. The shrinkage parameter is gradually annealed down over
time, and cross-validation (early stopping) is used for model
selection (regularization).

In this paper, MCHMMs have been shown to converge to local maxima in
likelihood space for a large non-parametric class of probability
density functions.  Empirical results, carried out in an artificial
and a gesture recognition domain, demonstrate that our approach
generalizes well even when trained with extremely scarce data.
Additional experiments characterize the natural trade-off between
sample set size and accuracy, illustrating that good results may be
achieved even from extremely small sample sets.  Additional results
described in \cite{Thrun98nABBREV} illustrate the effect of {\em
smoothing} in the forward/backward project step.

The choice of density trees for representing densities in MCHMMs is
somewhat ad hoc.  In principle, a range of methods could be applied,
such as Gaussian kernels. The trade-offs involved in choosing the right
density model are currently poorly understood and therefore warrant
future research.

% In all our experiments, MCHMM was trained to generate the most likely
% model, and discriminant analysis (Bayes rule) was applied for
% classification. If discrimination is the final task, it is unclear that
% this approach yields the best generalization; other training
% algorithms are therefore warranted.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






\section*{Acknowledgement}
\footnotesize
The authors gratefully acknowledge inspiring discussions with Frank
Dellaert, Nir Friedman, Daphne Koller, and Larry
Wasserman.

This research is sponsored in part by NSF (CAREER Award IIS-9876136)
and DARPA via TACOM (contract number DAAE07-98-C-L032), and Rome Labs
(contract number F30602-98-2-0137), which is gratefully acknowledged.
The views and conclusions contained in this document are those of the
authors and should not be interpreted as necessarily representing
official policies or endorsements, either expressed or implied, of
NSF, DARPA, TACOM, Rome Labs, or the United States Government.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\nocite{Gordon93aABBREV,Ghahramani98aABBREV,Briegel98aABBREV}




%\bibliography{../bib/my}\bibliographystyle{plain}
\begin{thebibliography}{10}
\footnotesize\parskip -1pt
\bibitem{Baum66aABBREV}
L.E. Baum and T.~Petrie.
\newblock Statistical inference for probabilistic functions of finite state
  markov chains.
\newblock {\em Annals of Mathematical Statistics}, 37, 1966.

\bibitem{Briegel98aABBREV}
T.~Briegel and V.~Tresp.
\newblock A Monte Carlo generalized EM-type algorithm for state and parameter
  estimation in nonlinear state space models.
\newblock  {\em NIPS-98}.

\bibitem{Charniak93aABBREV}
E.~Charniak.
\newblock {\em Statistical Language Learning}.
\newblock MIT Press, 1993.

\bibitem{Dean88aABBREV}
T.~L. Dean and M.~Boddy.
\newblock An analysis of time-dependent planning.
\newblock {\em AAAI-92}.

\bibitem{Dellaert99aABBREV}
F.~Dellaert, D.~Fox, W.~Burgard, and S.~Thrun.
\newblock Monte Carlo Localization For Mobile Robots
\newblock {\em ICRA-99}.

\bibitem{Dempster77aABBREV}
A.P. Dempster, A.N. Laird, and D.B. Rubin.
\newblock Maximum likelihood from incomplete data via the em algorithm.
\newblock {\em Journal of the Royal Statistical Society, Series B},
  39(1), 1977.

\bibitem{Fox99aABBREV}
D.~Fox, W.~Burgard, F.~Dellaert, and S.~Thrun.
\newblock Monte Carlo Localization: Efficient position estimation for mobile
  robots.
\newblock {\em AAAI-99}.

\bibitem{Ghahramani98aABBREV}
Z.~Ghahramani and S.~Roweis.
\newblock Learning nonlinear stochastic dynamics using the generalized EM
  algorithm.
\newblock  {\em NIPS-98}.

\bibitem{Gordon93aABBREV}
N.~Gordon.
\newblock Novel approach to nonlinear/non-gaussian bayesian state estimation.
\newblock {\em IEE Proceedings (F)}, 140(2), 1993.


\bibitem{Hannaford91aABBREV}
B.~Hannaford and P.~Lee.
\newblock Hidden markov model analysis of force torque information in
  telemanipulation.
\newblock {\em International Journal of Robotics Research}, 10(5),
  1991.

\bibitem{Henderson96aABBREV}
J.~Henderson, S.~Salzberg, and K.~Fasman.
\newblock Finding genes in human {DNA} with a hidden {M}arkov model.
\newblock {\em ISMB-96}.

\bibitem{Isard98aABBREV}
M.~Isard and A.~Blake.
\newblock Condensation: conditional density propagation for visual tracking.
\newblock {\em International Journal of Computer Vision}, 1998.

\bibitem{Juang85aABBREV}
B.H. Juang.
\newblock Maximum likelihood estimation for mixture multivariate stochastic
  observations of markov chains.
\newblock {\em AT\&T Technical Journal}, 64(6), 1985.

\bibitem{Kanazawa95aABBREV}
K.~Kanazawa, D.~Koller, and S.J. Russell.
\newblock Stochastic simulation algorithms for dynamic probabilistic networks.
\newblock {\em UAI-95}.

\bibitem{Koller98aABBREV}
D.~Koller and R.~Fratkina.
\newblock Using learning for approximation in stochastic processes.
\newblock {\em ICML-98}.

\bibitem{Krogh94aABBREV}
A.~Krogh, I.S. Mian, and D.~Haussler.
\newblock A hidden markov model that finds genes in e. coli dna.
\newblock {\em Nucleic Acids Research}, 22, 1994.


\bibitem{Liporace82aABBREV}
L.A. Liporace.
\newblock Maximum likelihood estimation for multivariate observations of markov
  sources.
\newblock {\em IEEE Transactions in Information Theory}, 28(5), 1982.

\bibitem{Liu98aABBREV}
J.~Liu and R.~Chen.
\newblock Sequential monte carlo methods for dynamic systems.
\newblock {\em Journal of the American Statistical Association}, 93, 1998.

\bibitem{McCallum98aABBREV}
A.~McCallum, R.~Rosenfeld, T.~Mitchell, and A.Y. Ng.
\newblock Improving text classification by shrinkage in a hierarchy of classes.
\newblock {\em ICML-98}.

\bibitem{McLachlan97aABBREV} G.J. McLachlan and T.~Krishnan.
\newblock {\em The {EM} Algorithm and Extensions}.  \newblock Wiley,
1997.

\bibitem{Meila96aABBREV}
M.~Meil\u{a} and M.I. Jordan.
\newblock Learning fine motion by markov mixture of experts.
\newblock {\em NIPS-96}.


\bibitem{Moore97aABBREV}
A.W. Moore, J.~Schneider, and K.~Deng.
\newblock Efficient locally weighted polynomial regression predictions.
\newblock {\em ICML-97}.

\bibitem{Neal93aABBREV}
R.M. Neal.
\newblock Probabilistic inference using markov chain monte carlo methods.
\newblock TR CRG-TR-93-1, University of Toronto, 1993.

\bibitem{Neal98aABBREV}
R.M. Neal and G.E. Hinton.
\newblock A view of the {EM} algorithm that justifies incremental, sparse, and
  other variants.
\newblock In {\em Learning in Graphical Models}. Kluwer, 1998.

\bibitem{Omohundro91aABBREV}
S.~M. Omohundro.
\newblock Bumptrees for Efficient Function, Constraint, and Classification
  Learning.
\newblock {\em NIPS-91}.


\bibitem{Pitt99aABBREV}
M.~Pitt and N.~Shephard.
\newblock Filtering via simulation: auxiliary particle filter.
\newblock {\em Journal of the American Statistical Association}, 1999.


\bibitem{Rabiner89aABBREV}
L.~R. Rabiner.
\newblock A tutorial on hidden markov models and selected applications in
  speech recognition.
\newblock In {\em Proceedings of the IEEE}, 1989.


\bibitem{Rubin88aABBREV}
D.B. Rubin.
\newblock Using the {SIR} algorithm to simulate posterior distributions.
\newblock In {\em Bayesian Statistics 3}. Oxford University Press, 1988.

\bibitem{Shatkay97aABBREV}
H~Shatkay and L.~Kaelbling.
\newblock Learning topological maps with weak local odometric information.
\newblock {\em IJCAI-97}.

\bibitem{Stein55aABBREV}
C.~Stein.
\newblock Inadmissibility of the usual estimator for the mean of a multivariate
  normal distribution.
\newblock In {\em Proceedings of the Third Berkeley Symposium on Mathematical
  Statistics and Probability 1}, 1995.

\bibitem{Tanner93aABBREV}
M.A. Tanner.
\newblock {\em Tools for Statistical Inference}.
\newblock 2nd ed., Springer, 1993.


\bibitem{Thrun98nABBREV}
S.~Thrun and J.~Langford.
\newblock Monte carlo hidden markov models.
\newblock TR CMU-CS-98-179, Carnegie Mellon University, 1997.

\bibitem{Waibel90aABBREV}
A.~Waibel and K.-F. Lee, editors.
\newblock {\em Readings in Speech Recognition}.
\newblock Morgan Kaufmann, 1990.

\bibitem{Zilberstein95aABBREV}
S.~Zilberstein and S.~Russell.
\newblock Approximate reasoning using anytime algorithms.
\newblock In  {\em Imprecise and Approximate Computation}.
  Kluwer, 1995.

\end{thebibliography}

\end{document}


 

