\documentstyle[12pt,times,ml99sub,fleqn,../styles/smallcaptions,my,epsf]{article}



\title{Sample-Based Hidden Markov Models\thanks{This paper has only been submitted to ICML-99}}
\author{Sebastian Thrun\ \ \  John C.~Langford\ \ \ Dieter Fox}
\date{}

\begin{document}

\maketitle

\begin{center}
School of Computer Science\\
Carnegie Mellon University\\
5000 Forbes Ave.\\
Pittsburgh, PA 15213-3891, USA
\end{center}

\begin{abstract}
\noindent
We present a learning algorithm for 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, just as in regular HMM
learning.  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.  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,
which illustrate the appropriateness of the approach in practice.
\end{abstract}

\noindent
{\bf Keywords:} {\em Hidden Markov Models, Sampling, Monte Carlo Methods, Shrinkage, Density Trees, Gesture Recognition}

\noindent
{\bf Email address of contact author:}  thrun@cs.cmu.edu

\noindent
{\bf Phone number of contact author:} +1 (412) 268-8077  

\noindent
{\bf Electronic version:} gzipped PostScript emailed to icml99@ijs.si
\newpage


\section{Introduction}




Throughout the last decade, hidden Markov models (HMMs) have been
applied successfully to a large range of temporal signal processing
problems, for tasks such as modeling, filtering, classification and
prediction of time sequences partially observable, 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 hidden Markov models with {\em 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 to choose the right
model.

MCHMMs differ from conventional HMMs in four primary aspects:%\\[-8mm]
\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. Conventional HMM can only support
discrete state and observation spaces (see
\cite{Liporace82aABBREV,Juang85aABBREV} for extensions to real-valued
observation spaces with parametric models). In reality, however,
processes modeled by HMMs are often continuous in nature.  For
example, the state space of the vocal tract of human beings, which
plays a primary role in the generation of speech, is continuous; yet
HMMs in speech recognition are usually discrete, even though they
attempt to model the speech-generating process
\cite{Rabiner89aABBREV,Waibel90aABBREV}.
\item
{\bf Not confined to specific parametric density.}  Virtually all
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, we typically lack such strong a
prior knowledge. 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,  MCHMMs can be used in an {\em any-time}
fashion. Any-time algorithm 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.%\\[-7mm]
% \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 describes MCHMMs. We present the MCHMM convergence theorem,
which proves with probability one that the Baum-Welch learning
algorithm converges 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 validate the approach in practice.





\section{Continuous Space 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 distributions are continuous, Lipschitz, non-zero over
some compact interval, and possess continuous density functions.



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, distribution
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
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}
\alpha_t^{(n)}(x)&=& \int \alpha_{t-1}^{(n)}(x') \;\mu^{(n)}(x\mid x')\; \nu^{(n)}(O_t \mid x) \;dx'
\label{eqAlpha}\\
\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') &=& \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}' }
\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)&&
\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{eqPi}
\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 decreases 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 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. %\hfill$\Box$



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




\section{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%\\[-7mm]
\begin{eqnarray}
X&=&\{\langle x_{1},p_{x_1}\rangle\, \ldots,\langle x_{N},p_{x_N}\rangle \}
%\;\;\;\;\mbox{with}\;\sum_{n=1}^{N}
p_{x_n}=1
\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 random variables over the event space
$\{x_{1},\ldots,x_{N}\}$, whose distributions are defined by
$\{p_{x_1}, \ldots, p_{x_N}\}$ \cite{Neal93aABBREV,Pitt99aABBREV}.
Sampling can be applied to approximate arbitrary probability densities
$f$ \cite{Tanner93aABBREV}.




%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.0mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,65)(0,0)
\put(1,17){\special{psfile=pictures/hmm-21.ps vscale=23 hscale=23}}
\put(1,-17){\special{psfile=pictures/hmm-20.ps vscale=23 hscale=23}}
\put(40,-40){\special{psfile=pictures2/hmm.1.pdf.rot.ps vscale=38 hscale=38}}
\put(3,63){\footnotesize (a)}
\put(3,29){\footnotesize (b)}
\put(50,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}
\ \\[-13mm]
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%

We will distinguish two methods for sampling:%\\[-7mm]
\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 $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 a convolved densities and other densities conditioned on
observations (c.f., (\ref{eqAlpha})). Instead of sampling directly
from $f$, importance sampling samples from a different distribution
$g$ that has the property that $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}
p_{x_n}&=& f(x_n)/g(x_n)
\end{eqnarray}
Importance sampling \cite{Rubin88aABBREV} has been studied extensively
in the statistics literature (see for example \cite{Tanner93aABBREV}).%\\[-7mm]
\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}).  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 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{Density Trees}
\label{SectDensTrees}

While sample sets are sufficient to approximate continuous-valued
distributions, they differ from those in that they are {\em
discrete}. This is problematic if one wants to {\em multiply}
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.

MCHMM solves this problem by transforming sample sets into
(non-pa\-ra\-me\-tric) 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:%\\[-7mm]
\begin{enumerate}\itemsep 0mm\parskip 0mm \parindent 0mm
\item 
At least $\sqrt{N}$ samples $\langle x,p_x\rangle\in X$ fall into in
$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 interval $v$ is divided into two equally sized
intervals 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}
Then 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).}



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

\section{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
estimations 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. If $\rho=0$, on the other hand,
there is no shrinking 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-------------------------------------------------------------------------%
\begin{table}[t]
\begin{scriptsize}
\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
 Copy the sample set representing $\pi$ into $\alpha_0$ (c.f.,
 (\ref{alphaInit})).
\item
Set $\beta_{T} = 1$.
\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{eqBeta})). For each $t$ with
  $1\leq t < T$ 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 $\langle x,p_x \rangle $ a
  probability 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 sampled $\langle x$ a probability $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.
\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.
\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:] Increase $N$.
\end{description}
\end{scriptsize}
\ \\[-18mm]
%\rule{\textwidth}{.2mm} 
\caption{The MCHMM algorithm at-a-glance.}
\label{TabMCHMM}
\ \\[-13mm]
\end{table}
%#@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}). A MCHMM is a computational instantiation of a GHMM that
represents all densities through samples and trees.  It applies
likelihood-weighted sampling for forward and backward projection
(c.f., Equations (\ref{eqAlpha}) and (\ref{eqBeta})), and it uses
annealing and cross-validation to determine the best shrinkage factor.
To ensure convergence, the number of samples $N$ is increased over time.

The learning algorithm for MCHMM is depicted in
Table~\ref{TabMCHMM}.  MCHMMs use both sample-based and tree
representations during learning. After learning, it suffices to store
only the tree-based version of the model $\lambda= \{\pi,\mu,\nu\}$;
all sample sets can be discarded. When applying a trained MCHMM to a
new data set (e.g., for analysis, prediction, or classification), only
the ``forward'' densities $\alpha_t(x)$ have to be estimated (just
like in conventional HMMs, Kalman filters, or dynamic belief networks
\cite{Dean89aABBREV}). For that, no trees have to be grown. Instead,
samples for $\alpha_{t+1}(x)$ are obtained by sampling from the sample
set representing $\alpha_{t}(x)$ (see also Section~\ref{SectSamp1})
and the tree representing $\mu$.  This is equivalent to the
``condensation algorithm'' which has been successfully applied in
domains such as computer vision and robotics
\cite{Dellaert99aABBREV,Isard98aABBREV}.

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 a MCHMM to converge to a local maximum
in with probability 1. }

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





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






\section{Experimental Results}
\label{Sect7}
We have applied MCHMMs to two problems, an artificial one
which was chosen for demonstration purposes, and a more difficult
real-world gesture recognition problem. 

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.
\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:
\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.
\end{eqnarray}      
other independent noise variables $\varepsilon$.


%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{0.8mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,36)(25,33)
\put(-6,-20){\special{psfile=pictures/gest-20.ps vscale=40 hscale=40}}
\put(36,-20){\special{psfile=pictures/gest-21.ps vscale=40 hscale=40}}
\put(78,-20){\special{psfile=pictures/gest-22.ps vscale=40 hscale=40}}
\put(120,-20){\special{psfile=pictures/gest-23.ps vscale=40 hscale=40}}
\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}
\ \\[-13mm]
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%


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.


%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.0mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,43)(0,0)
\put(-18,-41){\special{psfile=pictures/hmm-1c.ps vscale=45 hscale=45}}
\put(60,-41){\special{psfile=pictures/hmm-1.ps vscale=45 hscale=45}}
\put(1,40){\footnotesize (a)}
\put(79,40){\footnotesize (b)}
\end{picture}
\caption{Log score of 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
score of ``positive'' training data (same orientation as used for
training), whereas the bottom curve shows the log score 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-------------------------------------------------------------------------%




%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.0mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,43)(0,0)
\put(-18,-41){\special{psfile=pictures/hmm-3.ps vscale=45 hscale=45}}
\put(60,-41){\special{psfile=pictures/hmm-4.ps vscale=45 hscale=45}}
\put(1,40){\footnotesize (a)}
\put(79,40){\footnotesize (b)}
\end{picture}
\caption{(a) Log score for the gesture data base obtained for
independent testing data. The top curve shows the log score of
gestures drawn in same same way as the training gesture, whereas the
bottom curve shows the log score of gestures drawn in opposite
direction. (b) Trade-off between sample set size and log score, for
the gesture data set. The more samples are available, the higher the
score, and hence the more accurate the results. However, even
extremely small sample sets yield reasonable discrimination.}
%
\label{FigRes2}%key%
%\rule{\textwidth}{.2mm}
\ \\[-13mm]
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%




Figures~\ref{FigRes0} shows results obtained for the first dataset.
Shown in both figures are curves that characterize the ``log score''
as a function of the number of iterations. This score is the
real-valued analogue to the likelihood; however, since densities can
be larger than one, the score can also be larger than one, and hence
its logarithm is not bounded above.  Figure~\ref{FigRes0}a illustrates
that the log score 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 score for set of random sequence that are generated from the
same model as the training sequences.  The bottom curve depicts the
log score 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 score for
gestures drawn in the same direction as the training gesture, whereas
the bottom curve shows the log scores for gestures drawn in 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 score 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 score 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 score of gestures drawn in
the same direction as the training data, whereas the bottom curve
shows the log score of gestures drawn in opposite direction.
The result in Figure~\ref{FigRes2}b illustrates that the score (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 score space (hence it grows
in score). 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. The processing of such sample sets is several
orders of magnitude faster than the hand motion involved in gesture
generation.  This makes MCHMMs extremely fast for on-line tracking and
discrimination.


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



\section{Conclusion}
\label{Sect9}
We have presented a new algorithm for hidden Markov models, called
Monte Carlo Hidden Markov Models (MCHMM). MCHMMs extend HMMs to
real-valued state and observation spaces. 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 over-fit) 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 to select the best model complexity.  The pervasive
use of Monte Carlo sampling led to the design of an any-time
implementation, capable of dynamically trading off computational
complexity and the accuracy of the results. Consequently, MCHMMs are
well-suited for time-critical applications that are frequently
encountered in the real world.

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 domain and
a more challenging gesture recognition domain, demonstrate the 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.

We conjecture that MCHMMs are better suited for many real-world
application domains such as speech and robotics applications than
conventional HMMs, for primarily three reasons: their support of
continuous representations of observations and state, their built-in
mechanisms for model selection, which reduces the burden of picking
the ``right'' model (e.g., number of states in conventional HMMs),
and, finally, their support of any-time computation, which makes them
extremely compliant in time-critical applications with bounded
computational resources.


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





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

This research is sponsored in part by DARPA via AFMSC
(contract number F04701-97-C-0022), TACOM (contract number
DAAE07-98-C-L032), and Rome Labs (contract number F30602-98-2-0137).
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
DARPA, AFMSC, TACOM, Rome Labs, or the United States Government.



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






%\bibliography{../bib/my}\bibliographystyle{plain}
%\scriptsize\parskip -4pt
\begin{thebibliography}{10}
\footnotesize\parskip -2pt
\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{Charniak93aABBREV}
E.~Charniak.
\newblock {\em Statistical Language Learning}.
\newblock MIT Press, 1993.

\bibitem{Dean89aABBREV}
T.~Dean and K.~Kanazawa.
\newblock A model for reasoning about persistence and causation.
\newblock {\em Computational Intelligence}, 5(3), 1989.

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


 

