\documentstyle[fleqn,jair,times,twoside]{article}
\firstpageno{1}

\jairheading{??}{??}{\pageref{firstpage}--\pageref{lastpage}}{12/98}{??}

\ShortHeadings{Monte Carlo Hidden Markov Models}{S.~Thrun and J.~Langford}

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



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




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



\begin{document}
\label{firstpage}

\title{Monte Carlo Hidden Markov Models}
\author{Sebastian Thrun \and John Langford\\
       \addr{School of Computer Science, Carnegie Mellon University,
	Pittsburgh, PA 15213}}



% For research notes, remove the comment character in the line below.
% \researchnote

\maketitle


\begin{abstract}
We present a learning algorithm for hidden Markov models
with continuous state and observation spaces.  All necessary
probability density functions 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 obtained
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}






% \keywords{annealing, any-time algorithms, Baum-Welch, density trees,
% early stopping, EM, hidden Markov models, machine learning, maximum
% likelihood estimation, Monte Carlo methods, temporal signal
% processing}




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

\section{Introduction}


Over the last decade or so, hidden Markov models have enjoyed an
enormous practical success in a large range of temporal signal
processing domains.  Hidden Markov models are often the method of
choice in areas such as speech recognition
\cite{Lee90aABBREV,Lee89aABBREV,Rabiner89aABBREV}, natural
language processing \cite{Charniak93aABBREV}, robotics
\cite{Meila96aABBREV,Koenig96aABBREV,Shatkay97aABBREV}, biological
sequence analysis
\cite{Henderson96aABBREV,Krogh94aABBREV,Pedersen96aABBREV}, and time
series analysis \cite{Hannaford91aABBREV,Weigend93aABBREV}. They are
well-suited for modeling, filtering, classification and prediction of
time sequences in a range of partially observable, stochastic
environments.

With few exceptions, existing HMM algorithms assume that both the
state space of the environment and its observation space are
discrete. Some researchers have developed algorithms that support more
compact feature-based state representations
\cite{Ghahramani96aABBREV,Saul96aABBREV} which are nevertheless
discrete; others have successfully proposed HMM models that can cope
with real-valued observation spaces
\cite{Liporace82aABBREV,Juang85aABBREV,Shatkay97aABBREV}. Kalman
filters \cite{Kalman60aABBREV,Welch95aABBREV} can be thought of as
HMMs with continuous state and action spaces, where both the state
transition and the observation densities are linear-Gaussian
functions. Kalman filters assume that the uncertainty in the state
estimation is always {\em normally} distributed (and hence unimodal),
which is too restrictive for many practical application domains (see
e.g., \cite{Burgard98cABBREV,Isard98aABBREV}).

In contrast, most ``natural'' state spaces and observation spaces are
continuous.  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 trained to model the speech-generating process
are typically discrete.  Robots, to name a second example, always
operate in continuous spaces; hence their state spaces are usually
best modeled by continuous state spaces.  Many popular sensors
(cameras, microphones, range finders) generate real-valued
measurements, which are better modeled using continuous observation
spaces. In practice, however, real-valued observation spaces are
usually truncated into discrete ones to accommodate the limitations of
conventional HMMs. A popular approach along these lines is to learn a
{\em code-book} (vector quantizer), which clusters real-valued
observations into finitely many bins, and thus maps real-valued sensor
measurements into a discrete space of manageable size
\cite{Waibel90aABBREV}.  The discreteness of HMMs is in stark contrast
to the continuous nature of many state and observation spaces.
%However, we lack sufficiently general
%algorithms that can learn {\em continuous} HMMs, with
%continuous state and observations spaces.

Existing HMM algorithms possess a second deficiency, which is
frequently addressed in the AI literature, but rarely in the
literature on HMMs: they do not provide mechanisms for adapting their
computational requirements to the available resources. This is
unproblematic in domains where computation can be carried out
off-line. However, trained HMMs are frequently employed in
time-critical domains, where meeting deadlines is essential. {\em
Any-time} algorithms \cite{Dean88aABBREV,Zilberstein95aABBREV} address
this issue.  Any-time algorithms can generate an answer at any time;
however, the quality of the solution increases with the time spent
computing it.  An any-time version of HMMs would enable them to adapt
their computational needs to what is available, thus providing maximum
flexibility and accuracy in time-critical domains. Marrying HMMs with
any-time computation is therefore a desirable goal.

This paper presents {\em Monte Carlo Hidden Markov Models
(MCHMMs)}. MCHMMs employs continuous state and observation spaces, and
once trained, they can be used in an any-time fashion. Our approach
employs Monte Carlo methods for approximating a large, non-parametric
class of density functions. To combine multiple densities (e.g., with
Bayes rule), it transforms sample sets into density trees.  Since
continuous state space are sufficiently rich to overfit any data set,
our approach uses shrinkage as 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.  We prove the convergence of MCHMMs
to a local maximum in likelihood space.  This theoretical results
justifies the use of MCHMM in a wide array of applications.  In
addition, empirical results are provided obtained for an artificial
domain and a gesture recognition domain, which illustrates the
robustness of the approach in practice.


The remainder of this paper is organized as
follows. Section~\ref{Sect2} establishes the basic terminology for
generalizing HMMs to continuous state and observation spaces. We will
then, in Section Section~\ref{Sect3}, discuss commonly used sampling
schemes and provide an algorithm for generating piece-wise constant
density trees from these samples. A key result in this section is a
proof of asymptotic consistency of both sample-based and tree-based
representations.  Section~\ref{Sect4} describes our approach for
complexity control (regularization) using shrinkage, followed by the
statement of the MCHMM algorithm in Section~\ref{Sect5}.
Section~\ref{Sect6} proves the MCHMM convergence theorem, which states
that under mild assumptions MCHMMs converge with high probability.
Empirical results are described in Section~\ref{Sect7}, which
specifically investigates MCHMM in the finite sample case. The
empirical results show that MCHMMs work well even in the
non-asymptotic case. Section~\ref{Sect7} also provides an experiment
that characterizes the relation between sample set size (computation)
and accuracy. Finally, related work is discussed in
Section~\ref{Sect8} and the paper is summarized in
Section~\ref{Sect9}.





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



\section{Generalized Hidden Markov Models}
\label{Sect2}
This section introduces {\em generalized hidden Markov models} (in
short: {\em GHMM}). GHMMs generalize conventional hidden Markov models
(HMMs) in that all spaces, state and observation, are continuous. Our
description closely follows that of Rabiner \cite{Rabiner86aABBREV},
with densities replacing finite probability distributions
throughout. Throughout this paper, we assume all event spaces and
random variables (e.g., state, observations) are measurable. We also
assume that unless otherwise specified, all probability distributions
are continuous and possess continuous density functions. Further
below, when introducing density trees, we will also assume that
densities are non-zero over a compact, bounded region, and that they
obey a Lipschitz condition.

A GHMM is a partially observable, time-invariant Markov chain with
continuous state and observation spaces and discrete time.  Let $x$
denote a {\em state variable} (a measurable random variable), defined
over some continuous space (e.g., $\Re^k$ for some $k$). 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}.  Densities are measurable function over the set
of Borel sets, hence the Riemann integral
\begin{eqnarray}
\int_{x_0}^{x_1} \mu(x'\mid x) \;dx'
\end{eqnarray}
measures, for $x_0< x_1$, the probability $Pr(x_0\leq x' < x_1 \mid  x)$
that the state preceding  $x$ lies in $[x_0,x_1)$.

In HMMs (and thus in GHMMs), state cannot be observed. 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)$ . If $b_0<b_1$,
\begin{eqnarray}
\int_{b_0}^{b_1} \nu(b\mid x) \;db
\end{eqnarray}
measures the probability that the observation $b$ is in $[b_0, b_1)$,
given that the state of the HMM is $x$. Thus, a generalized HMM is
uniquely defined trough three densities:
\begin{eqnarray}
\lambda&=& \{\pi,\mu,\nu\}.
\label{eqLambda}
\end{eqnarray}
Putting computational limitations aside for the moment, knowledge of
$\lambda$ is sufficient to tackle a variety of interesting practical
problems:
\begin{itemize}
\item
Computing distributions over states and observations at arbitrary
points in time,
\item
Generating representative example trajectories in state and
observation space,
\item
Determining the
likelihood of example trajectories under an HMM, and
\item
Classifying data sequences generated by mixtures of labeled HMMs.
\end{itemize}
Algorithms for these problems are described in detail in
\cite{Rabiner86aABBREV}; 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,Rabiner86aABBREV} 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)}$. As
shown for example in \cite{McLachlan97aABBREV,Neal98aABBREV}, both
steps increase the data likelihood $Pr(d\mid \lambda)$, or leave it
unchanged if, and only if, a local maximum in the likelihood function
has been reached.

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{Rabiner86aABBREV}, 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)
\label{alphaInit}\\
\beta_{T}^{(n)}(x) &=& 1
\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)
\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.  For finite-state
finite-observation HMMs, all densities and variables ($\mu$, $\nu$,
$\pi$, $\alpha$, $\beta$, $\gamma$, and $\xi$) can be implemented by
vectors, and the Baum-Welch algorithm has been proven to converge
\cite{Baum66aABBREV}.  Juang \cite{Juang85aABBREV} 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.  In
particular, the class of ellipsoidal symmetrical densities includes
Gaussians.  A GHMM can be viewed as the limit as the number of
Gaussians is allowed to increase to infinity, and then the number of
states is allowed to increase to infinity in Juang's analysis.  Since
Gaussians are fine bump functions any continuous \( \nu(O|x) \) can 
be the limit.  
\begin{eqnarray}
\lim _{K \rightarrow \infty} \sum _{k=1}^K c_{kx} \;b_{kx}(O)
\end{eqnarray}
where \( \sum_k c_{kx} = 1 \) and \( b_{kx}(O) \) is a Gaussian.

Juang's analysis contains the following forms of manipulation:
 \begin{eqnarray}
 \int_Y
\int_X f(x,y) \;dx\; dy &=& \int_X \int_Y f(x,y) \;dy\; dx 
\end{eqnarray}
 \begin{eqnarray}
 \nabla _Y
\int_X f(x,y) \;dx  &=& \int_X \nabla _Y f(x,y) \; dx 
\end{eqnarray}

where \( f(x,y) \) is some function proportional to a density.
Differentiability of all densities is a sufficient assumption 
for these statements to hold true.

{\bf Theorem~1. (GHMM Convergence Theorem)} {\sl If all of the input
distributions \( \lambda \) to a GHMM possess density functions that
are Lipschitz and differentiable in all variables, then for almost all
(measure 1) points the steps outlined above do not decrease the
likelihood of the observation sequence. They do not improve the
likelihood of the observation sequence at each iteration if, and only
if, the input distributions are at a critical point (local maximum, minima, or
saddle point).}

\underline{Proof.} As sketched above. (see also
\cite{Baum66aABBREV,Juang86aABBREV,McLachlan97aABBREV}, and see
\cite{Shatkay98aABBREV} for an extension to certain real-valued
spaces).   \hfill$\Box$



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



\section{Density Approximation}
\label{Sect3}
This section describes sample-based and tree-based methods for density
approximation. The notion of {\em asymptotic consistency} is
introduced and results, along with error bounds, are given for a
popular sampling algorithm and a tree method. Throughout this section,
we assume that all density functions are centered on a compact and
bounded region in $\Re^k$ (for an arbitrary $k$) and that they are
Lipschitz.


\subsection{Sampling}
\label{SectSamp1}
{\em Samples} are values drawn from the domain of a density function, where
each sample is annotated by a non-negative real value
\cite{Liu98aABBREV,Neal93aABBREV,Winkler95aABBREV}. Sample sets are
(finite) sets of samples, annotated by numerical probability values.
More specifically, let $f$ be a probability density function, and let
$N$ denote a positive number (the cardinality of a sample set).  Then
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 \}
\;\;\;\;\mbox{with}\;\sum_{n=1}^{N}
p_{x_n}=1
\end{eqnarray}
where $x_{n}\in {\rm dom}(f)$ and $p_{x_n}\in [0,1]$ for all $n$ with
$1\leq n \leq N$.  Sample sets can be thought of as discrete
distributions over the event space $\{x_{1},\ldots,x_{N}\}$ and
probability distribution defined by $\{p_{x_1}, \ldots, p_{x_N}\}$
\cite{Neal93aABBREV,Pitt99aABBREV}.

One popular sampling method is called {\em importance sampling}, which was
originally introduced by Rubin
\cite{Rubin88aABBREV}. Importance sampling approximates a density $f$
by drawing samples $x$ from a distribution with density $g$, and
assigning a weight $p_x$ proportional to $\frac{f(x)}{g(x)}$. Obviously,
the density $g$ must be non-zero over the support of $f$, i.e.,
\begin{eqnarray}
f(x)>0 &\Longrightarrow g(x)>0
\end{eqnarray}
We will distinguish two special cases:
\begin{enumerate}
\item
{\bf Likelihood-weighted sampling.}  If $g\mbox{$\equiv$}f$, we will refer to
the sampling algorithm as {\em likelihood-weighted sampling}.  Here we
draw samples according to $f$ and assigns equal probability to each
sample in the sample set
\cite{Kanazawa95aABBREV}:
\begin{eqnarray}
x_{n} && \mbox{is drawn according to $f$}\nonumber\\
p_{x_n} &=& N^{-1}
\end{eqnarray}
In many cases, likelihood-weighted sampling can easily be implemented
using {\em rejection sampling} \cite{Neal93aABBREV}.
\item
{\bf Uniform sampling.} If $g$ is uniform over (a superset of) the
support of $f$, we will call the sampling algorithm {\em uniform
sampling}. This sampling algorithm draws values randomly from the
support of $f$ (which is assumed to be bounded), and assigns to each
value $x$ a probability proportional to its density $f(x)$:
\begin{eqnarray}
x_{n} && \mbox{is drawn uniformly from ${\rm dom}(f)$}\nonumber\\
p_{x_n} &=& f(x_{n})\;\left[\sum_{i=1}^{N} f(x_{i})\right]^{-1}
\end{eqnarray}
\end{enumerate}
Uniform sampling covers the domain of a density uniformly regardless
of the nature of the density function.   Likelihood-weighted sampling
populates the space according to $f$, so that the density of sample is
proportional to the density $f$.  Likelihood-weighted sampling can be
said to ``waste'' fewer samples in low-likelihood regions of $f$
\cite{Pitt99aABBREV}.  In other words, if $X_{\rm u}$ and $X_{\rm lw}$
are sample sets generated from the same distribution using uniform
sampling, and likelihood-weighted sampling, respectively, the
following holds true in expectation:
\begin{eqnarray}
\sum_{\langle x,p_x\rangle\in X_{\rm u}} f(x) &\leq& 
\sum_{\langle x,p_x\rangle\in X_{\rm lw}} f(x)
\end{eqnarray}
were the equality holds in expectation if and only if $f$ is uniform
\cite{Tanner93aABBREV}.  Henceforth, we will focus our attention on
likelihood-weighted sampling whenever we have an explicit
representation of a probability distribution (and importance sampling
otherwise).

%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{0.8mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,65)(-10,0)
\put(-10,-29){\special{psfile=pictures/hmm-21.ps vscale=35 hscale=35}}
\put(74,-29){\special{psfile=pictures/hmm-20.ps vscale=35 hscale=35}}
\put(0,60){(a)}
\put(84,60){(b)}
\end{picture}
\caption{(a) Data set (b) partitioned by a density tree.}
%
\label{FigSamplesAndTree}%key%
\rule{\textwidth}{.2mm}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%


Figure~\ref{FigSamplesAndTree}a shows a sample set drawn by
likelihood-weighted sampling from a distribution that resembles the
shape of a sine wave in 2D. 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$.  In practice, likelihood-weighted sampling is
often given preference over uniform sampling---specifically if the
target density is known to only populate (with non-zero measure) a
small subspace of its domain.

Both sampling methods can equally be applied to {\em sample from a
sample set} (resampling).  Let $X$ be a sample set. Under uniform
sampling, a new sample set is generated by randomly drawing samples
$\langle x,p_x
\rangle$ from $X$ with a uniform distribution (regardless of the $p_x$
values).  The $p_x$-values in the new sample set are then
re-normalized so that they add up to 1. Under likelihood-weighted
sampling, samples are drawn from $X$ according to the (discrete)
probability distribution induced by their $p_x$-values. Each sample is
then assigned the same probability.  Sampling from sample sets plays
an important role in the Monte Carlo HMM algorithm described below.




\subsection{Asymptotic Consistency}
\label{SectAsymCons}

The central idea of sampling is to {\em represent} density functions
by finite sample sets, which are computationally advantageous. So what
does it mean for a sample set $X$ to ``represent'' a density $f$?
Obviously, sample sets represent discrete distributions and thus cannot
represent distributions that possess densities (hence are continuous).
However, 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\;\;\;\;\;\;\mbox{w.p. 1}
\label{eqConv}
\end{eqnarray}
Recall the  system of half-open intervals is an inducing system for
the sigma algebra over the $\Re^k$; hence, it suffices to show
convergence for those sets. 

A key observation is that both sampling algorithms are asymptotically
consistent---as are all importance sampler with $f(x)>0\Longrightarrow
g(x)>0$ \cite{Tanner93aABBREV}. This is formalized for the more important
likelihood-weighted sampling technique in the following theorem, which
also provides a probabilistic error bound for the finite case.

{\bf Theorem~2.} {\sl The likelihood-weighted sampling method is asymptotic
consistent. For finite $N$, $\varepsilon > 0$ and $x_0<x_1$, the
likelihood that the error between a density $f$ and a sample $X$
is larger than $\varepsilon$ can be bounded as follows:
\begin{eqnarray}
Pr\left(\;\left| \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\right|\;> \varepsilon\;\right) &\leq& 2\;e^{-2\varepsilon^2 N}
\label{eqTheorem2}
\end{eqnarray}}
Thus, according to the theorem the error between $f$ and its sample
decreases at the familiar rate $\frac{1}{\sqrt{N}}$.

\underline{Proof.} The convergence follows directly from the Central
Limit Theorem (see also
\cite{Doucet98aABBREV,Pitt99aABBREV,Tanner93aABBREV}). The bound
follows from the Hoeffding bound, under the observation that for any
half-open interval $[x_0,x_1)$ each sample can be viewed as a zero-one
``guess'' of the size of the area $\int_{x_0}^{x_1}
f(x)\;dx$.\hfill$\Box$

More generally, the convergence rate of importance sampling is in
$O(\frac{1}{\sqrt{N}})$ if $frac{f(x)}{g(x)} < \infty$. The variance
of the error depends on the ``mismatch'' between $f$ and $g$ (see
e.g., \cite{Tanner93aABBREV}, pg.~33).  The \( N \) will be reduced by
a constant factor of \( \max_{x} \frac{f(x)}{g(x)} < \infty \) for
appropriate \( g(x) \).




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





\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},
that is, even in the limit they assign non-zero likelihood to only a
countable number of points. This is problematic if one wants to {\em
combine} densities represented through different sample sets: For
example, let $f$ and $g$ be two density functions defined over the
same domain, and let $X$ be a sample of $f$ and $Y$ a sample of
$g$. Then with probability 1, none of the samples in $X$ and $Y$ are
identical, and thus it is {\em not} straightforward how to obtain an
approximation of their product $f\cdot g$ from $X$ and $Y$. Notice
that multiplications of densities are required by the Baum-Welch
algorithm (see e.g., Equation (\ref{eqGamma})).

Density trees, which are quite common in the statistical literature
\cite{Koller98aABBREV,Moore97aABBREV,Omohundro87aABBREV,Omohundro91aABBREV},
transform sample sets into density functions.  Unfortunately, not all
tree growing methods are asymptotically consistent when applied to
samples generated from a density $f$.  We will describe a simple
algorithm which we will prove to be asymptotically consistent.

Our algorithm annotates each node in the tree with a hyper-rectangular
subspace of ${\rm dom}(f)$, denoted by $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$. A node $i$ is split whenever the
following two conditions are fulfilled:
\begin{enumerate}
\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$.
\end{enumerate}
If a node is split, its interval $v$ is divided into two equally sized
intervals along its longest dimension. These intervals are assigned to
the two children of the node. Otherwise, a node becomes a leaf node of
the density tree. 

%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.0mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,87)(0,0)
\put(0,-50){\special{psfile=pictures2/hmm.3.pdf.rot.ps vscale=50 hscale=50}}
\end{picture}
\caption{Density tree. Each branch cuts the area under a node into two
equally-sized, rectangular halves, shown in grey.}
%
\label{FigDensTree}%key%
\rule{\textwidth}{.2mm}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%

An example of such a density tree defined over $\Re^2$ is illustrated
in Figure~\ref{FigDensTree}. Here the root node covers the entire
domain of $f$. Each child covers exactly half the space of each parent
node. The areas covered by each leaf are mutually exclusive; their
union covers the entire space.  Figure~\ref{FigSamplesAndTree}b shows
a density tree generated for the sample set shown in
Figure~\ref{FigSamplesAndTree}a.  The reader may notice that the
criteria of the conditions above have been chosen to facilitate our
proof of asymptotic consistency; in practice, a much wider range of
criteria for stopping tree growth can be applied.


Density trees define density functions that continue discrete sample
sets into $\Re^k$. 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.

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

\subsection{Asymptotic Consistency of Density Trees}

An important result is the asymptotic consistency of density trees.
This result, which justifies the use of density trees are
non-parametric approximation, will be formalized in the
following theorem:

{\bf Theorem~3.} {\sl Density trees are asymptotically
consistent, that is, for all $x_0, x_1 \in{\rm dom(f)}$ with
$x_0<x_1$ and for any density tree $\hat{f}(x)$ generated from a
sample set $X$ generated from a density $f$, the following
holds true:
\begin{eqnarray}
\lim_{N \rightarrow \infty}\;\;
\int_{x_0}^{x_1} \hat{f}(x)\;dx
&= &
\int_{x_0}^{x_1} f(x)\;dx \;\;\;\;\;\;\mbox{w.p. 1}.
\label{eqTheorem3a}
\end{eqnarray}}
% If $N$ is finite, the following error bound applies
% \begin{eqnarray}
% Pr\left(\;\left|
% \int_{x_0}^{x_1} \hat{f}(x)\;dx
% - \int_{x_0}^{x_1} f(x)\;dx\right| > \varepsilon\right) &\leq & ...?
% \label{eqTheorem3a}
% \end{eqnarray}
% for arbitrary $\varepsilon > 0$.


We will  first prove a related result, namely that of
pointwise stochastic convergence
\begin{eqnarray}
\lim_{N \rightarrow \infty}\;\;
\hat{f}(x) &=
&f(x)\;\;\;\;\;\;\;\;\;\;(\forall x\in {\rm dom}(f))\;\;\;\;\;\;\mbox{w.p. 1}
\end{eqnarray}
which under the assumptions made in this paper implies the theorem.
The proof is carried out in two stages, each of which is described by
its own lemma.


{\bf Lemma~1.} {\sl The density $\hat{f}$ converges to a function
$\bar{f}$ which is defined through the same tree as $\hat{f}$, but
with the weighted relative frequencies $\hat{\sigma}_i$ replaced by
their true frequencies, denoted $\bar{\sigma_i}$:
\begin{eqnarray}
\bar{\sigma_i} &=&\int_{V_i} f(x)\;dx.
\end{eqnarray}}

{\bf Lemma~2.} {\sl The density $\hat{f}$ converges to $f$.}

\underline{Proof of Lemma~1.}  For each leaf $i$, the Hoeffding bound
states that
\begin{eqnarray}
Pr(|\bar{\sigma_i} - \hat{\sigma}_i| > \varepsilon )&<&
2\;e^{-2\varepsilon^2 N}.
\end{eqnarray}
Our tree growing rule limits the depth of the tree
to $\lfloor \frac{1}{4} \log_2 N \rfloor$. Hence, there are
at most 
\begin{eqnarray}
2^{\frac{1}{4} \log_2 N} &=&  N^{\frac{1}{4}}
\end{eqnarray}
nodes in the tree. Thus, the probability that there exists a leaf
whose empirical frequency $\hat{\sigma}_i$ deviates from the true
probability $\bar{\sigma_i}$ is bounded by
\begin{eqnarray}
Pr(\exists \;{\rm leaf} \;i: |\bar{\sigma_i} - \hat{\sigma}_i| > \varepsilon) &<&
2\;N^{\frac{1}{4}}\;e^{-2\varepsilon^2 N}.
\label{eq1}
\end{eqnarray}
The error of the empirical frequency estimates is now translated
into density errors. Recall that according to (\ref{eqTreeDens}),
the relative frequency is related to the density of a tree by
the volume $V_i$ of each leaf. Consequently,
\begin{eqnarray}
|\bar{f}(x)-\hat{f}(x)|
&=& |V_i|^{-1}\; | \bar{\sigma_i} - \hat{\sigma}_i|
\end{eqnarray}
Observing that the volume of the interval covered by a leaf $|V_i|$
is at least $N^{-\frac{1}{4}}$---which directly follows from the depth
limit imposed  when growing density trees---we obtain
\begin{eqnarray}
&&|\bar{f}(x)-\hat{f}(x)|
\;\;\leq\;\; N^{-\frac{1}{4}}\; | \bar{\sigma_i} - \hat{\sigma}_i|
\\
&\Longleftrightarrow&
N^{\frac{1}{4}}\;|\bar{f}(x)-\hat{f}(x)|
\;\;\leq\;\;  | \bar{\sigma_i} - \hat{\sigma}_i|
\end{eqnarray}
Substituting this into (\ref{eq1}) yields:
\begin{eqnarray}
Pr( |\bar{f}(x)-\hat{f}(x)| > N^{\frac{1}{4}}\varepsilon) &<& 2\;N^{\frac{1}{4}}\;e^{-2\varepsilon^2 N},
\end{eqnarray}
which with the substitution $\varepsilon' =
N^{\frac{1}{4}}\varepsilon$ leads to the following error bound between
$\hat{f}$ and $\bar{f}$:
\begin{eqnarray}
Pr( |\bar{f}(x)-\hat{f}(x)| > \varepsilon' )&<& 2\;N^{\frac{1}{4}}\;e^{-2\varepsilon'^2  \sqrt{N}}.
\label{eq2}
\end{eqnarray}
Obviously, the right hand-size of (\ref{eq2}) converges to 0 as
$N\rightarrow \infty$, which proves pointwise convergence of $\hat{f}
$to$ \bar{f}$. \hfill$\Box$

It remains to be shown that $\bar{f}$ converges of $f$.

\underline{Proof of Lemma~2.} It suffices to show that with high
likelihood, any leaf $i$ that covers an interval with non-zero
measure, i.e.,
\begin{eqnarray}
\bar{\sigma}_i& >& 0
\end{eqnarray}
will be split infinitely often as $N$ is increased.  The desired
convergence then follows directly from the fact that $f$ obeys a
Lipschitz condition. The proof is given for likelihood-weighted
sampling.

Let $i$ be a leaf node, and let ${\rm depth}(i)$ denote its depth.
Furthermore, let $n_i$ be the number of samples in node $i$:
\begin{eqnarray}
n_i &=&  \sum_{\langle x,p_{x}\rangle\in X}I_{x\in V_{i}} \;\;=\;\;
N\;\hat{\sigma}_i
\end{eqnarray}
The second equality exploits the assumption that samples are generated
by likelihood-weighted sampling.  Recall that the node $i$ is split if
$n_i\geq \sqrt{N}$, for sufficiently large $N$.  Without loss of
generality, let us assume that
\begin{eqnarray}
N &>& \max\left\{\frac{1}{(\bar{\sigma}_i-\varepsilon)^2}\;,\;2^{{\rm
depth}(i)} \right\}\;\;\;\;\;\;\;\;\mbox{and}\;\;\;\;\;\;\;\;\varepsilon <
\bar{\sigma}_i
\label{eqMax}
\end{eqnarray}
The second term in the max ensures that the depth limit $\lfloor
\frac{1}{4} \log_2 N \rfloor$ in the tree growing rule does not
restrict splitting node $i$. The other assumptions (\ref{eqMax}) imply
that
\begin{eqnarray}
\bar{\sigma}_i \;\; > \;\; \frac{1}{\sqrt{N}} + \varepsilon
\end{eqnarray}
The Hoeffding bound now yields the desired result:
\begin{eqnarray}
&&Pr( \bar{\sigma}_i - \hat{\sigma}_i  > \varepsilon) \;\;\leq\;\; e^{-2\varepsilon^2 N}\\
&\Longrightarrow&
Pr\left(\frac{1}{\sqrt{N}} + \varepsilon  - \hat{\sigma}_i  > \varepsilon\right) \;\;\leq\;\; e^{-2\varepsilon^2 N}\\
&\Longleftrightarrow&
Pr\left(\frac{1}{\sqrt{N}} > \hat{\sigma}_i\right) \;\;\leq\;\; e^{-2\varepsilon^2 N}\\
&\Longleftrightarrow&
Pr\left(\frac{1}{\sqrt{N}} > \frac{n_i}{N}\right) \;\;\leq\;\; e^{-2\varepsilon^2 N}\\
&\Longleftrightarrow&
Pr\left(\sqrt{N} > n_i\right) \;\;\leq\;\; e^{-2\varepsilon^2 N}
\end{eqnarray}
Thus, with high probability $n_i \geq \sqrt{N}$ and node $i$ is split.
The pointwise convergence of $\bar{f}$ to $f$ now follows from our
assumption that $f$ is Lipschitz.\hfill$\Box$


\underline{Proof of the Theorem~3.} We have already shown that
\begin{eqnarray}
\lim_{N \rightarrow \infty}\;\;
\hat{f}(x) &=
&f(x)\;\;\;\;\;\;\mbox{w.p. 1 and}\; \forall x\in {\rm dom}(f).
\end{eqnarray}
Theorem~3 follows trivially from the triangle inequality, since ${\rm
dom}(f)$ is bounded and Lipschitz (which implies that $f$ is bounded).
\hfill$\Box$

The termination conditions for growing trees (see itemized list in
Section~\ref{SectDensTrees}) 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}
\label{Sect4}
We will now resume our consideration of GHMMs.  The continuous nature
of the state space in GHMMs, if represented by arbitrary sample sets
and/or trees, can easily 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 samples in the data set. 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 (see also
\cite{Vapnik82aABBREV}). Shrinkage trees were introduced in
\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.


%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.0mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,90)(0,0)
\put(0,-50){\special{psfile=pictures2/hmm.1.pdf.rot.ps vscale=50 hscale=50}}
\end{picture}
\caption{Shrinkage for complexity control. The density of each node is
an exponentially weighted mixture of densities of the nodes along the
path to the leaf. Shown here is an example for node $m$, whose path
leads through nodes $i$, $j$, $k$, and $l$. The $\rho$-terms
describe the mixture coefficients for this example. If $\rho$ is 0,
only the leaf estimate is used. For larger values of $\rho$,
estimates from data beyond the leaf data are used to estimate the
target density.}
%
\label{FigShrinkage}%key%
\rule{\textwidth}{.2mm}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%

Figure~\ref{FigShrinkage} 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
\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
0). Thus, $\rho$ starts at is maximum value (at which nothing can
be learned, since every density is uniform), which 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]
\caption{The MCHMM algorithm at-a-glance.}
\label{TabMCHMM}
\begin{small}
\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}
\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}
\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}
\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}
\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}

\item[M-step:]
\ \\[-4mm]
\begin{enumerate}
\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{small}
\rule{\textwidth}{.2mm} 
\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 described in detail 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
\cite{Kalman60aABBREV,Maybeck79aABBREV}, or dynamic belief networks
\cite{Dean89aABBREV,Russell95aABBREV}). 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$. The
$p_x$-values of $\alpha_{t+1}(x)$ are determined using the tree
representing $\nu$. This recursive resampling technique, known as
sampling/importance resampling \cite{Rubin88aABBREV,Smith92aABBREV}
applied to time-invariant Markov chains, converges at the rate
$1/\sqrt{N}$ (if $T$ is finite, c.f., Section~\ref{SectAsymCons}).
Sampling/importance resampling, which will further be discussed in
Section~\ref{Sect8}, has been successfully applied in domains such as
computer vision and robotics \cite{Dellaert98dABBREV,Isard98aABBREV}.




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


\section{Error Analysis and Convergence Results}
\label{Sect6}

This section presents the major convergence result for MCHMM. Under
mild assumptions, MCHMMs can be shown to converge with probability 1
to models $\lambda$ that locally maximize likelihood; just like
conventional HMMs. The proof builds on the well-known results for
discrete HMMs (see g.e.,
\cite{Baum66aABBREV,Juang85aABBREV,Juang86aABBREV,McLachlan97aABBREV,Shatkay98aABBREV}),
and shows that if the sample size $N$ is sufficiently large, the
deviation between the MCHMMM and the corresponding GHMM can be bounded
arbitrarily tightly with high probability, exploiting the asymptotic
consistency of the two major approximations: samples and density
trees.


\subsection{Relative Error}

The process of Monte Carlo simulation introduces error into the
Baum-Welch algorithm.  In order to understand how this affects
convergence, we need to understand how an initial error will propagate
through an EM step.  This analysis will also cover errors introduced
in the middle of an EM step as long as such errors converge to 0 with
large samples, because errors introduced in mid-step will be
indistinguishable from an error introduced at the beginning of the
calculation.  

What we would like to prove is that a small absolute initial error
will imply a small absolute final error.  Unfortunately, this is not
possible because several calculations are normalized integrals of
products.  Consider the integral \( \int _{x} f(x) g(x) dx \).  If \(
f(x) = step(x - 1/2) \) and \( g(x) = 1 - f(x) \) then \( \int _{x}
f(x) g(x) \;dx = 0 \).  However, if \( \bar{f}(x) = f(x) + \varepsilon
\) and \( \bar{g}(x) = g(x) + \varepsilon \) are integrated, we get a
non-zero integral: \( \int _{x=0} ^1 \bar{f}(x) \bar{g}(x) dx =
\varepsilon \).  First, notice that the calculation of \( \alpha _t
(x) \) is of this form.  Assume, for the moment, that \( \beta _t (x)
= 1 \).  Then \( \gamma _t (x) = \alpha _t (x) / \int \alpha_t (x')
dx' \).  Consequently, if \( \alpha _t (x) \) was small (\( <
\varepsilon \)) everywhere, \( \bar{\gamma} _t (x) \) can potentially
be very far from \( \gamma _t (x) \). Furthermore, any allowable error
can produce an arbitrarily large error in the output.

However, if we have a small {\em relative} error, 
\begin{eqnarray}
 \frac{\bar{f}(x)
- f(x)}{f(x)} &<& \varepsilon 
\end{eqnarray}
 we will be able to prove that a small
initial error implies a small final error.  The restriction to small
relative error is significant because density trees only guarantee a
small absolute error.  An extra assumption must be placed on the
initial distribution in order to guarantee that the density tree
produced by sampling from the distribution will have small relative
error.

To simplify the proof that a small initial error produces a small
final error, we only consider transforming functions to first order,
\( f(\bar{x}) \sim f'(x) (\bar{x} - x) \).  There is no zeroth order
error because the limit as the relative error approaches \( 0 \) is \(
f(x) \), \( \lim _{\bar{x} \rightarrow x} f( \bar{x} ) = f(x) \).
Higher order error becomes small quickly for all transforming
functions which we consider.  Consequently, we will be able to state
that for all \( \varepsilon, \delta \), there exists some number of
samples, \( N \), s.t.  \( f( \bar{x}) - f(x) \leq f'(x)( \bar{x} - x
) + \varepsilon \) with probability \( 1 - \delta \).

Our analysis does not consider the shrinkage factor $\rho$, which, since
$\lim_{n\rightarrow \infty}\rho^{(n)}=0$, has asymptotically no
effect.

\subsection{Convergence of MCHMMs}

We will now state the central result of this section: the MCHMM
convergence theorem and an important corollary.  The proof of the
theorem will be developed throughout the remainder of this section.

{\bf Theorem~4 (MCHMM Convergence Theorem).} {\sl If approximation of
the underlying distributions by density trees causes only a finite
relative error which decreases to \( 0 \) as \( N\rightarrow \infty \)
and an EM step in a GHMM starting with \( \pi ^{(n)} ,\mu ^{(n)} , \nu
^{(n)} \) produces output distributions \( \pi ^{(n+1)}, \mu ^{(n+1)},
\nu ^{(n+1)} \) then for all \( \varepsilon , \delta \) there exists
an \( N \) s.t.  the output of an MCHMM iteration taking \( \pi ^{(n)}
,\mu ^{(n)} , \nu ^{(n)} \) as inputs with probability \( 1 - \delta
\) satisfies \( |\pi ^{(n+1)}(x)-\pi ^{'(n+1)}(x)|< \varepsilon \), \(
|\mu ^{(n+1)}(x'|x)-\mu ^{'(n+1)}(x'|x)|< \varepsilon \), and \( |\nu
^{(n+1)}(O_{t}|x)-\nu ^{'(n+1)}(O_{t}|x)|< \varepsilon \) where \( \pi
^{'(n+1)} \), \( \mu ^{'(n+1)} \), and \( \nu ^{'(n+1)} \) are the
output distributions of the MCHMM iteration. }

The proof of Theorem~4 will be presented below, after introducing a
collection of useful lemmas. However, Theorem~4 implies the following
important corollary:

{\bf Corollary.} {\sl Under the same assumptions as in Theorem~4, any
strictly monotonically increasing schedule of \( N \)'s will cause a
MCHMM to converge to a local maximum in with probability 1. }

\underline{Proof of the corollary.}  The convergence of GHMMs is a
straightforward extension of the well-known convergence proof for HMMs
as outlined in the proof of Theorem~1.
%  We spare the reader the
%derivation, and refer him instead to the rich literature on this topic
%\cite{Baum66aABBREV,Juang85aABBREV,Juang86aABBREV,McLachlan97aABBREV};
%see \cite{Shatkay98aABBREV} for an extension to certain real-valued
%spaces. 
The action of a GHMM iteration causes movement through the
space of distributions.  We can view this action of a GHMM iteration as
inducing a vector field in the (infinite dimensional) space of
distributions with a magnitude proportional to the change in
likelihood and a direction given by \( \lambda ^{(n)} (x,x') - \lambda
^{(n+1)} (x,x') \).  The action of an MCHMM with probability \( 1 -
\delta \) will have a result in a perturbed vector with the error
bounded by \( \varepsilon \).  The repeated action of a MCHMM with N
increasing will create a Cauchy sequence where the \( \varepsilon \)
decreases to \( 0 \) and \( \delta \) decreases to \( 0 \).  This
Cauchy sequence will limit to local maxima with probability 1.  \hfill$\Box$

It is important to notice that this corollary does {\em not} state
that an MCHMM and GHMM starting with the same distribution will
converge to the {\em same} local maximum.  Such a result will generally
not hold true, as the finiteness of $N$ might influence the specific
local maximum to which a MCHMM converges.

It is also interesting to note that the noise introduced by MCHMMs in
the convergence step may be beneficial.  The argument is informal
because the noise of an MCHMM has both systematic (due to
representation) and random (due to sampling) effects.  First, consider
a GHMM which accidentally starts out at a local minimum (or a saddle
point).  The GHMM might be stuck here, but a MCHMM will, with
probability 1, step off of the local minima and move towards a local
maxima.  In addition, consider a landscape containing many local
maxima on it's flanks.  A GHMM could easily become trapped in a local
maxima while an MCHMM with a low $N$ will be kicked out of the local
maxima with high probability.  In this way, an MCHMM could gain some
of the benefits of simulated annealing.


\subsection{Propagation of Error}

Before we can prove the MCHMM convergence theorem, we will need to
develop an analysis of the propagation of errors through transforming
functions.

Let us start by assuming that we have some initial error \(
\overline{\pi }(x)=\pi (x)+\Delta _{\pi }(x) \), \(
\overline{\mu}(x|x')=\mu (x|x')+\Delta _{\mu }(x,x') \) and \(
\overline{\nu}(O_{t}|x)=\nu (O_{t}|x)+\Delta _{\nu }(O_{t},x) \).
These errors can, for example, be the errors introduced by
approximating a distribution with density trees.  It will also be
convenient to talk about the relative error
\begin{eqnarray}
 \Delta ^{r}_{\pi }(x)&=&\left|\frac{\Delta _{\pi }(x)}{\pi (x)}\right| \\
 \Delta ^{r}_{\mu }(x|x')&=&\left|\frac{\Delta _{\mu }(x|x')}{\mu (x|x')}\right| \\
 \Delta ^{r}_{\nu }(O_{t}|x)&=&\left|\frac{\Delta _{\nu }(O_{t}|x)}{\nu (O_{t}|x)}\right| \\
 \Delta ^{r}_{f}(y)&=&\left|\frac{\Delta _{f}(y)}{f(y)}\right| 
\end{eqnarray}
and the largest relative error over the input range
\begin{eqnarray}
\Delta ^{r}_{\pi }&=&\max _{x}\Delta ^{r}_{\pi }(x) \\
\Delta ^{r}_{\mu }&=&\max _{x,x'}\Delta ^{r}_{\mu }(x|x') \\
\Delta ^{r}_{\nu }&=&\max _{x,O_{t}}\Delta ^{r}_{\nu }(O_{t}|x) \\
\Delta ^{r}_{f}&=&\max _{y}\Delta ^{r}_{f}(y) 
\end{eqnarray}
The propagation of errors through individual calculations of
Baum-Welch will be characterized by taking a derivative.  Let \(
\bar{x} = x + \Delta(x) \).  The definition of a derivative is \(
\frac{d}{dx}f(x)=\lim _{ \bar{x} \rightarrow x}\frac{f( \bar{x}
)-f(x)}{ \bar{x} -x} \).  As the difference between 2 distributions,
\( \bar{x} \) and \( x \) approaches \( 0 \), the difference between
the resulting calculations, \( f(\bar{x}) \) and \( f(x) \), will
approach \( f'(x) \Delta (x) \) implying \( \Delta _{f} (x) =
f(\bar{x}) - f(x) = f'(x) \Delta (x) \) in the limit that \( \Delta
(x)\rightarrow 0 \).  In particular,
\begin{eqnarray}
\Delta _{f+g}(x) &=& \Delta _{f}(x)+\Delta _{g}(x) 
\label{Property1}\\
 \Delta _{\frac{g}{f}}(x) &=& \frac{-g(x)\Delta _{f}(x)+f(x)\Delta _{g}(x)}{f(x)^{2}} 
\label{Property2}\\ 
\Delta _{fg}(x) &=& f(x)\Delta _{g}(x) + g(x)\Delta _{f}(x) 
\label{Property3}\\ 
\Delta _{\int f}&=&\int \Delta _{f} 
\label{Property4} 
\end{eqnarray}
The last statement is only true for integration over compact regions
with finite valued \( f(x) \), because \( \int  \bar{f}(x) - f(x) dx
= \int \bar{f}(x) dx- \int f(x) dx\) under these assumptions.  These
assumptions were already made to prove convergence of density trees so
they are not an additional restriction.



\subsection{Error Bounds for Auxiliary Variables}
Given the above rules and notation we can calculate how the error will
propagate through the calculation of a \( \alpha _{t}(x) \).  The
following lemmas establish bounds on the relative errors of all
quantities calculated in a Baum-Welch iteration given some initial
error.  These lemmas are necessary to prove the MCHMM convergence
theorem.

{\bf Lemma~3.}
{\sl \( \Delta _{\alpha _{t}}^{r}\leq (t-1)(\Delta ^{r}_{\nu }+\Delta _{\mu }^{r})+\Delta ^{r}_{\pi } \) to first order. }

\underline{Proof.}
According to the property (\ref{Property4}),
\begin{eqnarray}
\Delta _{\alpha _{t}}&=&\Delta _{\int \alpha _{t-1}\mu \nu }(x)=\int _{x} \Delta _{\alpha _{t-1}\mu \nu }(x) \;dx
\end{eqnarray}
Using (\ref{Property3}), this expression is equivalent to
\begin{eqnarray} 
&=&\int _{x} \mu (x)\nu (x)\Delta _{\alpha _{t-1}}(x)+\nu (x)\alpha _{t-1}(x)\Delta _{\mu }(x)+\mu (x)\alpha _{t-1}(x)\Delta _{\nu }(x) \;dx
\nonumber\\
&=&
\int _{x}\mu (x)\nu (x)\alpha _{t-1}(x)\left(\frac{\Delta _{\alpha _{t-1}}(x)}{\alpha _{t-1}(x)}+\frac{\Delta _{\nu }(x)}{\nu (x)}+\frac{\Delta _{\mu }(x)}{\mu (x)}\right) \;dx 
\end{eqnarray}
Notice that the integral can be bounded by a maximum,
\begin{eqnarray} 
\Delta _{\alpha _{t}}&\leq& \alpha _{t}(x)\max _{x}\left(\left|\frac{\Delta _{\alpha _{t-1}}(x)}{\alpha _{t-1}(x)}+\frac{\Delta _{\nu }(x)}{\nu (x)}+\frac{\Delta _{\mu }(x)}{\mu (x)}\right|\right) 
\end{eqnarray} 
and,
\begin{eqnarray} 
\Delta _{\alpha _{t}}&\geq& -\alpha _{t}(x)\max _{x}\left(\left|\frac{\Delta _{\alpha _{t-1}}(x)}{\alpha _{t-1}(x)}+\frac{\Delta _{\nu }(x)}{\nu (x)}+\frac{\Delta _{\mu }(x)}{\mu (x)}\right|\right) 
\end{eqnarray} 
Division  by \( \alpha _{t}(x) \)and application
of the triangle inequality  yields the relative error bounds
\begin{eqnarray} 
\frac{\Delta _{\alpha _{t}}}{\alpha _{t}(x)} &\leq& \max _{x}\left(\left|\frac{\Delta _{\alpha _{t-1}}(x)}{\alpha _{t-1}(x)}\right|+\left|\frac{\Delta _{\nu }(x)}{\nu (x)}\right|+\left|\frac{\Delta _{\mu }(x)}{\mu (x)}\right|\right)
\\
\frac{\Delta _{\alpha _{t}}}{\alpha _{t}(x)} &\geq& - \max _{x}\left(\left|\frac{\Delta _{\alpha _{t-1}}(x)}{\alpha _{t-1}(x)}\right|+\left|\frac{\Delta _{\nu }(x)}{\nu (x)}\right|+\left|\frac{\Delta _{\mu }(x)}{\mu (x)}\right|\right)
\end{eqnarray} 
Consequently, we have that \( \Delta _{\alpha _{t}}^{r}\leq \Delta _{\alpha _{t-1}}^{r}+\Delta _{\nu }^{r}+\Delta _{\mu }^{r} \) 
  The equation is recursive, terminating for \( t = 0 \) with 
\( \Delta _{\alpha _{1}}^r = \Delta _{\pi}^r \).  At each step in the recursion, 
at most \(\Delta _{\nu }^{r}+\Delta _{\mu }^{r} \) is added to the error which implies 
\[
\Delta _{\alpha _{t}}^{r}\leq (t-1)(\Delta ^{r}_{\nu }+\Delta _{\mu }^{r})+\Delta ^{r}_{\pi }\]
which proves Lemma~3.
\hfill$\Box$

The calculation of the relative \( \beta _{t} \) error is analogous,
and the resulting bound it stated in the following lemma:

{\bf Lemma~4.}
{\sl \( \Delta _{\beta _{t}}^{r}\leq (T-t)(\Delta _{\nu }^{r}+\Delta _{\mu }^{r}) \) to first order. }

The proof is omitted, since it is similar to the proof for \( \Delta
_{\alpha _{t}}^{r} \). 
% Notice that Lemma~4 assumes a perfect
% representation for the uniform distribution.  An imperfect
% representation will add an extra (similarly bounded) error term to the
% bound.

The calculation of the error of \( \gamma _{t} \) is more difficult.

{\bf Lemma~5.}
{\sl \( \Delta _{\gamma }^{r}\leq 2(\Delta ^{r}_{\alpha }+\Delta ^{r}_{\beta })=2(\Delta ^{r}_{\pi }+(T-1)(\Delta ^{r}_{\mu }+\Delta ^{r}_{\nu })) \) to first order. }

\underline{Proof.} To simplify the notation, we will assume a \( t \)
subscript in the following proof. The results are independent of \( t
\).
\begin{eqnarray}
\Delta _{\gamma }(x)&=&\Delta _{\frac{\alpha \beta }{\int \alpha \beta }}(x)
\nonumber\\
&=&\frac{\Delta _{\alpha \beta }(x)\int _{x} \alpha (x)\beta (x) dx 
- \alpha (x)\beta (x)\Delta _{\int \alpha \beta }(x)}{\left(\int _{x} \alpha (x)\beta (x)\right)^{2}}  \nonumber\\
\nonumber\\
&=&\frac{[\beta (x)\Delta _{\alpha }(x)+\alpha (x)\Delta _{\beta }(x)]\int _{x} \alpha (x)\beta (x) dx}{(\int _{x} \alpha (x)\beta (x) dx)^{2}}\nonumber\\
&&
\hskip 15mm
-\frac{\alpha (x)\beta (x)\int _{x} [\alpha (x)\Delta _{\beta }(x)+\beta (x)\Delta _{\alpha }(x)]dx }
{(\int _{x} \alpha (x)\beta (x) dx)^{2}}\nonumber\\
\nonumber\\
&\leq&
\frac{[\beta (x)\Delta _{\alpha }(x)+\alpha (x)\Delta _{\beta }(x)]\int _{x} \alpha (x)\beta (x)dx}{(\int _{x} \alpha (x)\beta (x)dx)^{2}}\nonumber\\
&&
\hskip 15mm
+\frac{\alpha (x)\beta (x)[\int _{x} \alpha (x)\beta (x) dx]\max _{x}\left(\left|\frac{\Delta _{\beta }(x)}{\beta (x)}\right|+\left|\frac{\Delta _{\alpha }(x)}{\alpha (x)}\right|\right)}{(\int _{x} \alpha (x)\beta (x)dx)^{2}}\nonumber\\
\nonumber\\
&\leq&
\alpha (x)\beta (x)\frac{\left[\frac{\Delta _{\alpha }(x)}{\alpha (x)}+\frac{\Delta _{\beta }(x)}{\beta (x)}\right]+\max _{x}\left(\left|\frac{\Delta _{\beta }(x)}{\beta (x)}\right|+\left|\frac{\Delta _{\alpha }(x)}{\alpha (x)}\right|\right)}{\int _{x}\alpha (x)\beta (x) dx}
\end{eqnarray}

Through similar manipulations, we get:
\begin{eqnarray}
\Delta _{\gamma }(x)
&\geq& - 
\alpha (x)\beta (x)\frac{\left[\frac{\Delta _{\alpha }(x)}{\alpha (x)}+\frac{\Delta _{\beta }(x)}{\beta (x)}\right]+\max _{x}\left(\left|\frac{\Delta _{\beta }(x)}{\beta (x)}\right|+\left|\frac{\Delta _{\alpha }(x)}{\alpha (x)}\right|\right)}{\int _{x}\alpha (x)\beta (x) dx}
\end{eqnarray}
Now, we can divide by a factor of \( \gamma \) to get
\begin{eqnarray}
\Delta _{\gamma }^{r}&\leq& 2(\Delta ^{r}_{\alpha }+\Delta ^{r}_{\beta }) 
\end{eqnarray}
and apply Lemmas 3 and 4 for any value of \( t \) to get
\begin{eqnarray}
\Delta _{\gamma }^r &\leq& 2(\Delta ^{r}_{\pi }+(T-1)(\Delta ^{r}_{\mu }+\Delta ^{r}_{\nu }))
\end{eqnarray}
which completes the proof.
\hfill$\Box$


The bounds for \( \Delta ^{r}_{\xi } \) are similar to those of \(
\Delta _{\gamma }^{r} \), as documented by the following lemma.

{\bf Lemma~6.}
{\sl \( \Delta ^{r}_{\xi }\leq  \)\( 2(\Delta _{\pi }^{r}+(T-1)(\Delta _{\mu }^{r}+\Delta _{\nu }^{r})) \) to first order. }

The proof of Lemma~6 is analogous to that of Lemma~5 and therefore
omitted.  Notice that the bounds of \( \Delta _{\xi }^{r} \) and \(
\Delta _{\gamma }^{r} \) are independent of \( t \).

The M step calculations are done similarly. Let \( \pi ' \) = the next \( \pi  \) distribution,
\( \mu ' \) =the next \( \mu  \) distribution, and \( \nu ' \) = the next \( \nu  \) distribution. Then the following lemmas provide bounds
for the error of \( \pi' \), \( \mu' \), and \( \nu' \).

{\bf Lemma~7.}
{\sl \( \Delta ^{r}_{\pi '}=2(\Delta _{\pi }^{r}+(T-1)(\Delta _{\mu }^{r}+\Delta _{\nu }^{r})) \) to first order. }

\underline{Proof.}
Lemma~7 follows from Lemma~5, since
\( \Delta ^{r}_{\pi '}=\Delta _{\gamma _{0}}^{r} \).

{\bf Lemma~8.}
{\sl \( \Delta _{\mu '}^{r}=4(\Delta _{\pi }^{r}+(T-1)(\Delta _{\mu }^{r}+\Delta _{\nu }^{r})) \) to first order. }

\underline{Proof.}  
By property (\ref{Property2}), 
\begin{eqnarray}
\Delta _{\mu '}(x'|x)
&=&\Delta _{\frac{\sum \xi }{\sum \gamma }}(x'|x)\nonumber\\
&=&\left[-\left(\sum _{t=1}^{T - 1}\gamma _{t}(x)\right)\Delta _{\sum \xi }(x,x')+\left(\sum _{t=1}^{T - 1}\xi _{t}(x,x')\right)\Delta _{\sum \gamma }(x)\right]\frac{1}{\left(\sum _{t=1}^{T - 1}\gamma _{t}(x)\right)^{2}} \nonumber\\
&=& \left[-\left(\sum _{t=1}^{T - 1}\gamma _{t}(x)\right)\left(\sum _{t=1}^{T - 1}\Delta _{\xi _{t}}(x,x')\right)
\right.
\nonumber\\
&&\hskip 30mm \left.+\left(\sum _{t=1}^{T - 1}\xi _{t}(x,x')\right)\left(\sum _{t=1}^{T - 1}\Delta _{\gamma _{t}}(x)\right)\right]\frac{1}{\left(\sum _{t=1}^{T - 1}\gamma _{t}(x)\right)^{2}} \nonumber\\
&=&\left[\sum _{t=1}^{T - 1}\xi _{t}(x, x')\right]\;\left[-\frac{\sum _{t=1}^{T - 1}\Delta _{\xi _{t}}(x,x')}{\sum _{t=1}^{T - 1}\xi _{t}(x,x')}+\frac{\sum _{t=1}^{T - 1}\Delta _{\gamma _{t}}(x)}{\sum _{t=1}^{T - 1}\gamma _{t}(x)}\right]\;\frac{1}{\sum _{t=1}^{T - 1}\gamma _{t}(x)} \nonumber\\
&\leq& \left[\sum _{t=1}^{T - 1}\xi _{t}\right]\;\max _{t}\left(\frac{\Delta _{\xi _{t}}}{\xi _{t}}+\frac{\Delta _{\gamma _{t}}}{\gamma _{t}}\right)\frac{1}{\sum _{t=1}^{T - 1}\gamma _{t}} 
\end{eqnarray}
and, analogously,
\begin{eqnarray}
\Delta _{\mu '}(x'|x)
& \geq& - \left[\sum _{t=1}^{T - 1}\xi _{t}\right]\;\max _{t}\left(\frac{\Delta _{\xi _{t}}}{\xi _{t}}+\frac{\Delta _{\gamma _{t}}}{\gamma _{t}}\right)\frac{1}{\sum _{t=1}^{T - 1}\gamma _{t}} 
\end{eqnarray}
These equations are of the form: 
\begin{eqnarray}
\Delta _{\mu '}(x'|x)
& \leq& \mu(x'|x) \max _{t}\;\left(\frac{\Delta _{\xi _{t}}}{\xi _{t}}+\frac{\Delta _{\gamma _{t}}}{\gamma _{t}}\right)
\end{eqnarray}
and
\begin{eqnarray}
\Delta _{\mu '}(x'|x)
& \geq& - \mu(x'|x) \max _{t}\;\left(\frac{\Delta _{\xi _{t}}}{\xi _{t}}+\frac{\Delta _{\gamma _{t}}}{\gamma _{t}}\right)
\end{eqnarray}

which implies that \( \Delta _{\mu '}^{r}\leq (\Delta ^{r}_{\xi }+\Delta ^{r}_{\gamma }) \). The lemma follows from the 
bounds for \( \Delta ^{r}_{\gamma } \) and \( \Delta ^{r}_{\xi } \)
stated in Lemma~5 and 6.
\hfill$\Box$

{\bf Lemma~9.}
{\sl \( \Delta _{\nu '}^{r}=4(\Delta _{\pi }^{r}+(T-1)(\Delta _{\mu }^{r}+\Delta _{\nu }^{r})) \) to first order. }

The proof is analogous to the proof of Lemma~8.

{\bf Theorem~5.}
{\sl The error introduced in a round of Baum-Welch with some initial \( \Delta _{\pi }^{r},\Delta _{\mu }^{r}, \) and \( \Delta _{\nu }^{r} \) to first
order is \( \Delta ^{r}_{\pi '}=2((T-1)(\Delta ^{r}_{\nu }+\Delta ^{r}_{\mu })+\Delta ^{r}_{\pi }) \), \( \Delta ^{r}_{\mu '}=4((T-1)(\Delta ^{r}_{\nu }+\Delta ^{r}_{\mu })+\Delta ^{r}_{\pi }) \), and \( \Delta ^{r}_{\nu '}=4((T-1)(\Delta ^{r}_{\nu }+\Delta ^{r}_{\mu })+\Delta ^{r}_{\pi }) \). }

\underline{Proof.} Theorem~5 follows directly from
Lemmas 7 to 9.
\hfill$\Box$





\subsection{Implications for MCHMMs}

The analysis of error propagation for GHMMs implies that as \(
N\rightarrow \infty \) and the relative errors approach \( 0 \) the
error of a step in the generalized Baum-Welch algorithm tends towards
the first order errors: \( \forall \varepsilon >0:\exists N: \)\( \Delta
^{r}_{\pi '}=2((T-1)(\Delta ^{r}_{\nu }+\Delta ^{r}_{\mu })+\Delta
^{r}_{\pi })+\varepsilon \), \( \Delta ^{r}_{\mu '}=4((T-1)(\Delta
^{r}_{\nu }+\Delta ^{r}_{\mu })+\Delta ^{r}_{\pi })+\varepsilon \), and
\( \Delta ^{r}_{\nu '}=4((T-1)(\Delta ^{r}_{\nu }+\Delta ^{r}_{\mu
})+\Delta ^{r}_{\pi })+\varepsilon \). These first order errors also
converge to \( 0 \).
\begin{eqnarray}
\lim_{N\rightarrow \infty }\Delta _{\nu '}^{r}&=&0\\
\lim_{N\rightarrow \infty }\Delta _{\mu '}^{r}&=&0\\
\lim_{N\rightarrow \infty }\Delta _{\pi '}^{r}&=&0
\end{eqnarray}
We are now ready to present the central result of this section,
a proof of convergence for MCHMMs.

\underline{Proof of Theorem~4 (MCHMM Convergence Theorem).}  According to
Theorem~5, the first order relative errors of \( \pi ^{'(n+1)} \), \(
\mu ^{'(n+1)} \), and \( \nu ^{'(n+1)} \) are linear in the relative
errors of \( \pi ^{'(n)} \), \( \mu ^{'(n)} \), and \( \nu ^{'(n)}
\). By assumption, in the limit as \( N\rightarrow \infty \) these
errors go to \( 0 \) with probability 1, implying the theorem.
\hfill$\Box$


%\subsection{Error Analysis Conclusion}

MCHMMs converge to a local maximum, under the additional assumption
that the relative error of all distributions converge to \( 0 \) as \(
N\rightarrow \infty \).  Notice that the \( \rho \) parameter was not
used here---the \( \rho \) parameter selects between various models,
while we only proved convergence within one model. However, since
$\rho \longrightarrow 0$, its effect vanishes over time.

The MCHMM Convergence Theorem establishes the soundness of our
algorithm, and shows that MCHMMs can indeed be applied for learning a
large range of non-parametric statistical models with real-valued
state and observation spaces.  Empirical results using finite sample
sets, described in the next section, suggest that the algorithm is
stable even for small sample set sizes, and indeed maximizes data
likelihood in a Monte Carlo-fashion.




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






\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 experiments address the
following questions:
\begin{itemize}
\item Do MCHMMs converge empirically? If so, how good are the
resulting MCHMMs?
\item What accuracies can be obtained when using MCHMMs for discrimination?
\item How does the sample set size affect computational and accuracy
trade-offs?
\end{itemize}


The first data set, called the {\em noisy oscillation dataset},
consists of multiple sequences that basically oscillate around two
points. Observations are 10-dimensional and governed by
\begin{eqnarray}
O_t &=&\left\{\begin{array}{ll}
(0.25+\varepsilon_{t,1}, 0.25+\varepsilon_{t,2},0.25+\varepsilon_{t,3},\ldots, 0.25+\varepsilon_{t,10}) &
\mbox{if $t$ odd} \\
(0.75+\varepsilon_{t,1},0.75+\varepsilon_{t,2},0.75+\varepsilon_{t,3}, \ldots, 0.75+\varepsilon_{t,10}) &
\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}
(0.25+\varepsilon_{t,1}, 0.75+\varepsilon_{t,2},0.25+\varepsilon_{t,3}, \ldots, 0.75+\varepsilon_{t,10}) &
\mbox{if $t$ odd} \\
(0.75+\varepsilon_{t,1},0.25+\varepsilon_{t,25},0.75+\varepsilon_{t,3} \ldots, 0.25+\varepsilon_{t,10}) &
\mbox{if $t$ even} 
\end{array}\right.
\end{eqnarray}      
using new, independent noise variables. Notice that despite the fact
that there is a good amount of noise in the data, these data sets are
relatively easy to discriminate, since their observations fall into
different regions in the $\Re^{10}$.


%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{1.0mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,36)(28,33)
\put(-8,-20){\special{psfile=pictures/gest-20.ps vscale=50 hscale=50}}
\put(30,-20){\special{psfile=pictures/gest-21.ps vscale=50 hscale=50}}
\put(68,-20){\special{psfile=pictures/gest-22.ps vscale=50 hscale=50}}
\put(106,-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}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%



The second data set consisted of a collection of hand-drawn gestures,
represented in $\Re^2$. 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}{0.8mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,85)(-13,0)
\put(-13,-73){\special{psfile=pictures/hmm-1c.ps vscale=65 hscale=65}}
\end{picture}
\caption{Log score of training data as a function of the iteration of
EM, for the synthetic data set. 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$.  These graphs illustrate that MCHMMs tend to maximize the data
likelihood, even in the finite case. Each result is averaged over 10
independent experiments; 95\% confidence bars are also shown}
%
\label{FigRes0}%key%
\rule{\textwidth}{.2mm}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%
%\bigskip
%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{0.8mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,85)(-13,0)
\put(-13,-73){\special{psfile=pictures/hmm-1.ps vscale=65 hscale=65}}
\end{picture}
\caption{Log score of testing data as a function of the iteration of
EM, for the synthetic data set. The top curve shows the log score of
data generated from the model used in training, whereas the bottom
curve shows the log score for data generated from a different
model. 95\% confidence bars are also shown.}
%
\label{FigRes1}%key%
\rule{\textwidth}{.2mm}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%


%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{0.8mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,85)(-13,0)
\put(-13,-73){\special{psfile=pictures/hmm-3.ps vscale=65 hscale=65}}
\end{picture}
\caption{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.}
%
\label{FigRes2}%key%
\rule{\textwidth}{.2mm}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%


%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{0.8mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,85)(-13,0)
\put(-13,-73){\special{psfile=pictures/hmm-4.ps vscale=65 hscale=65}}
\end{picture}
\caption{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{FigRes3}%key%
\rule{\textwidth}{.2mm}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%

%#@a-------------------------------------------------------------------------%
\setlength{\unitlength}{0.8mm}
\begin{figure}[t]
%\rule{\textwidth}{.2mm} 
\begin{picture}(100,85)(-13,0)
\put(-13,-73){\special{psfile=pictures/hmm-8.ps vscale=65 hscale=65}}
\end{picture}
\caption{MCHMMs with (black curve) and without (grey curve) density
tree smoothing of $\alpha$ and $\beta$, under otherwise equal
conditions.  Each curve plots the log score of the testing data as a
function of the number of Baum-Welch iterations.  The grey curve is
equivalent to the top curve in Figure~\ref{FigRes2}.  These results
illustrate the negative effect of smoothing on log score of the
testing data.}
%
\label{FigRes4}%key%
\rule{\textwidth}{.2mm}
\end{figure}
\setlength{\unitlength}{1.0mm}
%#@b-------------------------------------------------------------------------%



Figures~\ref{FigRes0} and \ref{FigRes1} show 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.

The curves in Figure~\ref{FigRes0} show average log score for the {\em
training set} as a function of the number of iterations, averaged over
10 independent experiments (each curve).  The different curves were
obtained using different numbers of samples which, contrary to the
theoretical results, were kept constant throughout the experiments:
\begin{center}
\begin{tabular}{l|c|c}
Sample set size $N$ for$\ldots$ & $\mu,\nu$ & $\alpha,\beta,\gamma,\pi$\\
\hline
top curve: & 1,000 & 1,000 \\
middle curve: & 1,000 & 100 \\
bottom curve: & 100 & 10
\end{tabular}
\end{center}
Figure~\ref{FigRes0} also provides 95\% confidence bars for these
results.  The difference in absolute levels to which the scores appear
to converge stem from the fact that the maximum tree depth, and hence
the maximum of the tree density $\hat{f}$, is a function of the sample
set size; hence, if the number of samples is small, the tree will
remain shallow. The key observation, however, is their monotonicity.
These curves demonstrate that each iteration in MCHMM indeed increases
the data likelihood (or leaves it unchanged), even if only finitely
many samples are used. This finding, which we also consistently
observed for other data sets and parameter settings, illustrate that
MCHMMs are well-behaved in practice, that is, the finiteness of the
sample set appears not to cause catastrophic problems.



Figure~\ref{FigRes1} shows results for independent testing data.  The
upper curve in Figure~\ref{FigRes1} 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.  $N=100$ samples were used throughout for all
densities $\alpha$, $\beta$, $\gamma$, and $\pi$, and $N=1,000$
samples were used for $\mu$ and $\nu$, to account for their higher
dimensionality. The initial shrinkage factor was $\rho=1$, which was
annealed down at the rate $\bar{\rho}=0.9$.  In all experiments, the
dimension of the state space was $k=2$.  The results obtained for
lower-dimensional observation spaces, which are not shown here, were
even better. In all cases, the score clearly reflected that the MCHMM
had learned a good model. The classification accuracy for test sets was
consistently 100\%.


The result in Figure~\ref{FigRes1} illustrate the discriminative power
of MCHMM for this data set. It also demonstrates the effect of
annealing: The testing data likelihood is maximal at iteration 6 (thus
the optimal $\rho$ is $\hat{\rho}=0.53$), beyond which the MCHMM
starts overfitting the data.  In our experiments, the optimal stopping
point was consistently found within $\pm 1$ step, using a single,
independently generated sequence for cross-validation. In all our
experiments, the MCHMM consistently discriminated the two different
classes without error, regardless of the settings of the individual
parameters (we did not run experiments with $N<10$). The bars in
Figure~\ref{FigRes1} are confidence bars at the 95\% confidence level.

Figure~\ref{FigRes2} shows the corresponding 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}, the difference between classes
is smaller than in Figure~\ref{FigRes1}---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}
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.

A key advantage of MCHMM over HMMs lies in the ease of trading off
computational complexity and accuracy. In conventional HMM, the
computational complexity is fixed for any given (trained) model. In
contrast, MCHMM permit variable sampling set sizes during run-time,
after the model is trained. This leads to a nice {\em any-time}
algorithm \cite{Dean88aABBREV,Zilberstein95aABBREV}. Once trained, a
MCHMM can control its computational complexity on-line by dynamically
adjusting the sample size, depending on the available resources. The
any-time property is important for many real-world applications, where
computation resources are bounded. To the best of our knowledge, the
any-time nature is a unique property of MCHMMs, which is not shared by
previous HMM algorithms.

Figure~\ref{FigRes3} illustrates the trade-off between computation and
accuracy empirically for the gesture data base.  Shown there is the
trade-off between the number of samples and the likelihood of the
testing data. 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}, the
top curve in Figure~\ref{FigRes2} 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{FigRes3} 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{FigRes3}
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 in this domain. 

Following suggestions by Koller\footnote{personal communication}, we
evaluated the utility of {\em smoothing} in the context of
MCHMMs \cite{Doucet98aABBREV,Koller98aABBREV}.  Smoothing
transforms those sample sets that represent posteriors over states
($\alpha$ and $\beta$) into more compact representations such as
density trees, which are then used for likelihood-weighted sampling
instead of the original samples.  To apply smoothing to MCHMMs, the
basic learning algorithm (see Table~\ref{TabMCHMM}) was modified as
follows:
\begin{description}
\item[{\rm Step 3a:}]
  Generate $N$ samples $\langle x',p_{x'}\rangle$ from the {\em tree}
  representing $\alpha_{t-1}$ using likelihood-weighted sampling.
\item[{\rm Step 4a:}]
  Generate $N$ samples $\langle x',p_{x'}\rangle$ from the {\em tree}
  representing $\beta_{t+1}$ using likelihood-weighted sampling.
\end{description}
Notice that this is a minor modification of the basic algorithm, as
trees are readily computed for the various densities $\alpha$ and
$\beta$; but now they are also used for likelihood-weighted sampling.
Obviously, smoothing reduces the variance of the various $\alpha$ and
$\beta$ at the cost of increasing bias.  It has been suggested that
smoothing with density trees can further improve performance
\cite{Koller98aABBREV}.

Unfortunately, in our experiments we were unable to confirm the
utility of smoothing in the context of MCHMMs. For the artificial data
set, the generalization performance of MCHMMs with smoothing was
indistinguishable from that obtained without smoothing, both in terms
of log score of the testing data and the class discrimination
accuracy. Smoothing actually worsened the performance for the gesture
data set. Figure~\ref{FigRes4} shows the log score of testing data
using smoothing (black curve), and compares it with the corresponding
log scores obtained without smoothing (gray curve, copied from
Figure~\ref{FigRes2}). Smoothing was found to reduce the log score of
the testing data; thus, MCHMMs trained with smoothing are less capable
of ``explaining'' new examples. This is not surprising, as smoothing
is not an information loss-free operation, and it is not clear that
the inductive bias provided by a tree representation is beneficial.
The reduction of the log score led to a slight reduction of the
discrimination accuracy. Using cross-validation, the average
classification accuracy for MCHMMs with smoothing was 82.0\%, which is
4.0\% smaller than the accuracy obtained without smoothing.



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

\nocite{Baldi94aABBREV}

\section{Related Work}
\label{Sect8}

Hidden Markov Models \cite{Baum66aABBREV,Rabiner86aABBREV} have been
successfully applied to a huge range of applications requiring
temporal signal processing.  As indicated in the introduction of this
paper, most state-of-the-art speech recognition systems rely on HMMs
(see, for example, the various papers in \cite{Waibel90aABBREV}).
Recently, extensions of HMMs have been successfully applied in the
context of mobile robotics
\cite{Koenig96aABBREV,Shatkay97aABBREV,Thrun98cABBREV}, where they
have led to improved solutions of the {\em concurrent mapping and
localization} problem, which is generally acknowledged as one of the
most difficult problems in robotics
\cite{Cox91aABBREV,Kortenkamp97aABBREV}.  These are
just two examples of successful application areas of HMM; the number
of successful applications is numerous.

Most HMM algorithms differ from MCHMMs in that they only apply to
discrete state and observation spaces. While the majority of HMM
approaches represent states individually (in a flat way), some more
recent approaches have extended HMMs to loosely coupled factorial
representations \cite{Ghahramani96aABBREV,Saul96aABBREV}, which
represent state by features but are nevertheless discrete. Others have
successfully proposed HMM models that can cope with real-valued
observation spaces \cite{Liporace82aABBREV,Juang85aABBREV}.  Our
approach generalizes HMMs to continuous state and observation spaces
with a large, non-parametric range of state transition and observation
densities. It also differs from previous HMM approaches in that it
provides a mechanisms for complexity control (regularization); in
previous approaches, the complexity of the internal state space had to
be calibrated carefully by hand. Finally, MCHMMs provide a mechanisms
to trade off computational requirements and accuracy in an any-time
fashion, which differs from previous approaches whose computational
requirements were not adjustable during run-time
\cite{Dean88aABBREV,Zilberstein95aABBREV}.  We envision that all of
these advances are important for a broad range of practical problems.

Sampling techniques, one of the two methods used for density
approximation in MCHMMs, have recently gained popularity in the
applied statistics and AI literature. Various researchers have applied
sampling methods in the context of state estimation in dynamical
systems
\cite{Boyen98aABBREV,Doucet98aABBREV,Kanazawa95aABBREV,Liu98aABBREV,Pitt99aABBREV}
and learning \cite{Freitas98aABBREV}. A nice introduction into the use
of sampling for density approximation can be found in
\cite{Neal93aABBREV}.  Our approach to state estimation (computation
of the $\alpha$s) is essentially equivalent to the {\em condensation
algorithm} proposed by Isard and Blake \cite{Isard98aABBREV} and the
{\em Markov localization algorithm} proposed by Dellaert and
colleagues \cite{Dellaert98dABBREV,Dellaert98eABBREV}, both of which
are basically versions of the well-known {\em sampling/importance
resampling} (SIR) algorithm
\cite{Rubin88aABBREV,Smith92aABBREV}. Similar approaches are known as
{\em particle filters} \cite{Pitt99aABBREV}, {\em bootstrap}
\cite{Efron93aABBREV}, and {\em survival of the fittest}
\cite{Kanazawa95aABBREV}. All these approaches are concerned with
state estimation in an HMM-like or Kalman filter-like fashion.  Thus,
they rely on the a priori availability of $\mu$ and $\nu$, which are
learned from data by the MCHMM algorithm proposed here.  To our
knowledge, the application of sampling-based methods to learning
non-parametric state transition models and observation models in HMMs
is new, as is our proposal for the integration of a forward ($\alpha$)
and backward ($\beta$) phase using trees. The theoretical results in
this paper demonstrate that our approach can be applied to a large
class of problems, assuming that sufficiently many samples are used.
Trees have frequently been used for density approximation, most
recently in
\cite{Koller98aABBREV,Moore97aABBREV,Omohundro87aABBREV,Omohundro91aABBREV}. However,
we are not aware of a proof of asymptotic consistency for density
trees, although we suspect that such a result exists.






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


\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. Both representations were proven to be
asymptotically consistent under minimal assumptions on the nature of
the density that is being approximated.  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.

We have proved the asymptotic consistency of MCHMMs for a large 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.

While the theoretical results of this paper hold only in the limit,
actual applications have to live with finite sample sets. Thus, an
interesting and open question is under what conditions MCHMM with
finite sample sets converges, and what error to expect. While
principal error bounds of this type are relatively easy to obtain in
static approximation setting, the recursive nature of Baum-Welch makes
it difficult to obtain bounds for MCHMM. The empirical results
described here, however, suggest that MCHMM is a robust method that
works well in practice even with relatively 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}
%\footnotesize\parskip 2pt



\label{lastpage}
\end{document}


 
