\documentclass[12pt]{article}
\usepackage{latex2e_dswf}
\usepackage{epsfig,graphicx}
\usepackage{varioref,chicago}

\newcommand{\LTE}{\mbox{\textsc{l}}_{\mbox{\scriptsize\textsc{2}}}%
  \mbox{\textsc{e}}}
\newcommand{\LTEs}{\mbox{\textsc{l}}_{\mbox{\scriptsize\textsc{2}}}%
  \mbox{\textsc{e}}}
\newcommand{\lte}{\hbox{\scshape {\itshape L$_2$E}}\,}
\newcommand{\bx}{{\bf x}}
\newcommand{\bmu}{{\mbox{\boldmath$\mu$}}}

\title{The Stochastic Mode Tree and Clustering}
\author{David W.~Scott and William F.~Szewczyk\footnote{David W.~Scott
is Professor, Department of Statistics, MS-138, Rice University,
Houston, TX 77005-1892 (Email:  scottdw\symbol{64}stat.rice.edu).
William F.~Szewczyk is a senior mathematician with the National Security
Agency (Email:  wfszewc\symbol{64}afterlife.ncsc.mil).}}
\date{\today}  %% January 6, 14, 2000

\begin{document}
\fancytitle

\begin{abstract}
\input{abs}
\end{abstract}
	
\section{Introduction}

The search for structure in an unknown set of data
$\{\bx_i\in\Re^d,\,i=1,\ldots,n\}$ is one of the most
challenging problems in statistics.  Graphical approaches
can reveal important details about individual variables or
pairs of variables, such as skewness, correlation, and clumping
\cite{xgobi}.  %% (XGobi, 199x).
However, no graphical exploration can be expected to
uncover all structure in higher dimensions.  Without any prior
information on a parametric density form, a nonparametric approach such
as the multivariate kernel estimator
\begin{equation}
{\hat f}(\bx)=\frac{1}{n\,|H|} \sum_{i=1}^n
	K \left( H^{-1}(\bx-\bx_i)\right)\,,
	\label{mkde}
\end{equation}
may be contemplated \cite{mvde}.  %%% (Scott, 1992).
Once the matrix of smoothing parameters $H$ is calibrated,
the kernel estimate can be probed
for the presence of simple forms of structure.  The most
important features on the density surface are modes.
Indeed, modes may naturally be associated with clusters
%% \cite{cox70,goodgaskins,hartigan,sager82}. %% (Sager, 198x).
\cite{hartigan75,goodgaskins,sager82}. %% (Sager, 198x).

In general, cluster analysis provides a sophisticated toolbox for the discovery
of structure in data, but clustering algorithms do not always rely
upon probability models.  Hierarchical clustering algorithms
bypass the density function to form clusters iteratively.
Starting with each data point as an individual cluster,
the ``closest'' pair is combined into a new cluster recursively.
A history of the process may be presented graphically in a binary
tree known as the {\it dendrogram}; for example, see Figure \ref{singlink:dendro}
for a dendrogram of the male blue crab data \cite{ripley}.
While almost entirely assumption- and model-free, the tree
is not unique, depending heavily upon the choice of metric,
definition of distance between two sets of points (clusters),
and any initial transformation of the data.
Finally, the dendrogram embodies many possible clusterings,
and the task of choosing where to cut the tree remains.
Without some model assumptions, this last task must remain
largely {\it ad hoc}.

\begin{figure}[htbp]
 \begin{center}
  \epsfig{file=FIGS/fig01.eps,width=4.0in,angle=270}
    \caption{
	Single linkage dendrogram of
	sphered/blurred male blue {\it Leptograpsus} crab data.}
    \label{singlink:dendro}
 \end{center}
\end{figure}

An appealing flexible class of probability models is the mixture
of Normals (MON),
\begin{equation}
  {\hat f}(\bx)=\sum_{k=1}^K w_k\,\phi(\bx|\bmu_k,\Sigma_k)\,,
  \label{def:mixnor}
\end{equation}
where $\phi(\cdot|\bmu,\Sigma)$ denotes a multivariate normal
density with mean $\bmu$ and covariance $\Sigma$.
When fitted to data by the EM algorithm \cite{dlr77}, %%% (Titterington, et al., 19xx),
this model provides a powerful and easily interpretable
summary of the structure in
the data by associating clusters with components.
Cluster centers are given by the means, $\bmu_k$,
cluster shapes by the covariance matrices, $\Sigma_k$,
and cluster sizes by the weights, $w_k$.
Observe the close relationship between
MON and the kernel method in Equation (\ref{mkde}), which is a highly
constrained MON (with all the data ``rounded'' to one of the $K$
values, $\bmu_k$, and all $\Sigma_k=H^2$).

In practice, a precise
choice of $K$ is not possible for high-dimensional data.
In theory, any multivariate density
may be represented by an infinite mixture of normals, but the
representation is by no means unique in the $L_2$ sense
(since the basis functions are not orthogonal).  One cluster
might be well-represented by a single multivariate normal, while another
cluster may require dozens of normal components to capture
skewness (and yet still be unimodal).  The non-uniqueness grows
exponentially with $d$ and is yet another manifestation of
the curse of dimensionality.

Furthermore, the stability of the estimates of the
parameters of the mixture decreases as $K$ increases
(although the overall density estimate may be quite stable).
Care must be exercised in the EM algorithm to avoid the
introduction of Dirac spikes into the final estimate.
An engineering solution is to take $K$ much larger than
necessary, carefully control the progress of the EM
algorithm, so as to be relatively confident in the overall
quality of $\hat f$ but with no expectation of
interpretability of the estimated parameters.
For example, one model of NIST speech data features $K=2048$
mixtures in $d=39$ dimensions \cite{reynolds95}.
Post-analysis \cite{fktm} %%% (Scott and Szewczyk, 2000)
showed that $K$ could be reduced by more than half with no change in
the performance of $\hat f$ for discrimination.
Of course, the parameters when $K=1024$ are hardly more interpretable!

\citeN{br93}  %% Banfield et al.~(199x)
have carefully constructed a hierarchy of mixture models
by allowing $\Sigma_k$ to be an identity matrix, $I_d$,
an arbitrary diagonal matrix, proportional to each other, or fully
unconstrained.  The authors use BIC \cite{bayes_factor} to determine $K$.
Such an approach generally results in a small $K$, but in practice has
trouble with large dimensions ($d>10$).  The paradox is that as $d$
increases, $\Sigma_k$ should be less constrained,
but practical issues require the $\Sigma_k$ to be {\it more} constrained.

The relationship between the
number of modes in a mixture density and $K$ is a complicated
issue; see \citeN{tsm}.  Obviously, the mixture density (\ref{def:mixnor})
may have many fewer than $K$
modes.  However, \citeN{scott.modes.97} %%% Scott/Szew (1997)
showed that the mixture (\ref{def:mixnor})
may have {\it more} than $K$ modes when $d>1$,
%% Thus the earlier strategy of probing the surface
%% of a kernel estimator for its modes may provide important
%% information when mixture models cannot be fit, or are overfit,
%% or are poorly fit.
but any extra modes tend to be very small.
\citeN{boswell83} %%% Boswell (198x)
successfully probed the kernel surface
in $d=100$ dimensions.  Numerical optimization can locate
almost all sample modes in the kernel estimate.
However, finding all the modes is expensive and no test exists that assures
every mode has been found.

With kernel estimates, the use of the same smoothing parameter
everywhere in $\Re^d$ leads to heavy oversmoothing in some
regions and heavy undersmoothing in other regions.
Indeed, \citeN{terrell.scott92} %%% Terrell and Scott (1987)
show how $H$ should vary locally (i.e., covariance shape vary as well as size).
Thus in practice, any fixed kernel estimate may be expected to suppress
modes in rough areas (oversmoothing) and to display false modes
near peaks or tails in smooth areas (undersmoothing).
\citeN{minnotte93}   %%% Minnotte and Scott (1993)
found a simple algorithm that avoids this local adaptivity
problem.  They suggest finding the locations of all modes as a function of the
smoothing parameter.  When $d=1$ and $K=\phi$,
the resulting figure, which they called a {\it mode tree},
looks remarkably like a dendrogram;  see Figure \ref{mt.1D:dendro}.
(Note that there are only 36 unique body diameter values among the 50 cases.)
The special properties of the Normal kernel guarantee that the
mode tree is a binary tree \cite{silverman81}.
When $h\approx0$ each unique data point is its own mode.  As $h$ increases, pairs of
modes collapse into each other until eventually a single final mode exists
near $x={\bar x}=13.35\hbox{mm}$.

\begin{figure}[htbp]
 \begin{center}
  \epsfig{file=FIGS/fig01.5.eps,width=4.0in,angle=270}
    \caption{
	Mode tree of the body diameter (mm) of 50 male
	blue {\it Leptograpsus} crabs.}
    \label{mt.1D:dendro}
 \end{center}
\end{figure}

When $d>1$, with $H=h\cdot I_d$, the authors show that the mode
tree may still be formed.
%% (modulo small inversions due to the
%% non-monotonicity of the number of modes as $h$ increases).
A smooth mode tree exists in $\Re^{d+1}$ but is difficult to
display when $d>2$.  A simple algorithm converts the mode tree
into a dendrogram (where the horizontal axis has no particular
relationship to the original variables), with the ``dissimilarity''
measured by the bandwidth, $h$.  For example,
the mode tree using all 5 variables of the 50 male blue crabs is displayed
in Figure \ref{mt.5D:dendro}.  The result is somewhat similar to
the single linkage clustering in Figure \ref{singlink:dendro}.

 While no single kernel estimate with fixed bandwidth
can be expected to show all of the real modes and only the
real modes, the family of kernel estimates is likely to contain
all of the modes.  %%% , except very small modes when $n$ is also small.
``Real'' modes tend to persist (vertically)
in the mode tree, while ``false'' modes cluster around real
modes and seldom persist very long before splitting.
Of course, singleton outliers ($x=20\hbox{mm}$, e.g.) are exceptions,
and the use of the logarithmic vertical scale greatly magnifies
the finer details.
\citeN{minnotte97} %%% Minnotte (1997, AOS)
demonstrated one way to ``test'' the veracity of a mode at each split;
see also \citeN{hartigan}, \citeN{sun97}, \citeN{marron98},
\citeN{scott.bumps.98}, and \citeN{walthers99}.  %%  and \citeN{ieee00}.
%%% (Other refs:  IEEE, Marron testing, J Sun, Walthers Stanford, 199x).
The mode tree approach captures the essence of a successful
density-based clustering strategy but is still limited in the testing arena
beyond one dimension.

\begin{figure}[htbp]
 \begin{center}
  \epsfig{file=FIGS/fig02.eps,width=4.0in,angle=270}
    \caption{
	Mode tree of the sphered/blurred male blue {\it Leptograpsus} crab data.}
    \label{mt.5D:dendro}
 \end{center}
\end{figure}

The purpose of this paper is to extend recent results
on {\it single-component} mixture fitting by a minimum-distance
criterion to a {\it stochastic mode tree} that allows probing of
high-dimensional data at multiple resolutions, without actually
constructing a complicated locally adaptive estimator \cite{scott.l2e.99}.
Given the present limitation that a fully parameterized locally adaptive estimator is
nearly impossible to calibrate with today's technology for $d>3$,
the approach outlined here is of fundamental interest.


\section{The Mode Tree}

The mode tree is a natural attempt to link probability
modeling with hierarchical clustering.
The construction of the mode tree from the bottom up is easy
to describe.  When $h\approx0$, the kernel estimate will have
a mode at each unique data point.  As the bandwidth $h$ is
sequentially increased ($1-2\%$ at each step), the change in the location
of each mode is estimated.  The resultant list of modes is examined
for ties (that is, two modes have collapsed) and that information
is recorded in the binary tree.  The process continues until
the final pair of modes collapses into a single mode.

In practice, care must be exercised.  If too large a
change in the bandwidth $h$ is used, then several modes may appear
to collapse at one step, and the exact linkages may be recorded
in error.  This may be avoided by decreasing the
change in $h$ until two and only two modes collapse.
But the surprisingly subtle aspect in this step can be the
difficulty of a successful search
for the new mode location as the bandwidth is increased.
Because one of the modes that merges first becomes a saddle
point, and then ``jumps'' to the appropriate adjacent mode,
Newton's method is likely to take too great an initial step and
converge to the ``wrong'' mode.  Thus the search should maintain a
local hill-climbing behavior.
Yet, strict local hill-climbing algorithms are well-known
to be very slow to converge.
In independent work, \citeN{wong.mt} %%% Wong (199x)
proposed a functional iteration algorithm to update modal locations.
While this algorithm appears parameter free, a careful
examination of its performance with data in two dimensions shows that the
algorithm does not always converge to the correct adjacent mode
(as determined by not wishing to cross gradient contours).
Typically such an error does not have dire consequences, as these modes
would have eventually joined to form a cluster and are
simply being combined in the wrong order further down the tree.  In practice, we may
hope that such subtle errors do not occur very often.

For many iterations, the small change in bandwidth hardly
changes the location of any modes. Much of the computational effort, however, is
spent dealing with such (essentially stationary) modes.
To determine which modes are not changing, one could calculate the
Hessian at each mode.  (All of the Hessian's eigenvalues must be
negative, of course, at a local maximum.)
Then one could simply observe which
modes have any eigenvalues that are approaching zero (indicating
that the mode is about to go ``saddle'' and jump to some
nearby mode).  Computation could then be focused only on those
nearly ``critical'' modes.  However, the amount of additional
computation %%% may exceed updating the modal locations.  Such effort
might be justified in difficult cases.

For large datasets, it would be more convenient to construct
the mode tree from the top down, since detail near the bottom of
the mode tree is not important.  However,
as the bandwidth is decreased, finding the location of
the newly emerging mode can
be difficult.  Typically, when a mode ``splits,'' one of the new modes
remains at the current location.  The other mode ``jumps'' somewhere
in the vicinity, but finding it among the dozens of other
modes nearby may or may not succeed, especially when $d>>1$.
Furthermore, if the increment on
the bandwidth $h$ is too great, then several new modes may have
appeared, and there is no way to know or prove that all new modes have been found.
Thus, while the bottom-up approach may be slow, its
behavior is much more understandable and controllable.

Once a ``correct'' mode tree is computed, the decision
about which modes are ``real'' and which are just noise begins.
\citeN{minnotte97} %% Minnotte (1997)
proposed a local bootstrap procedure to test the
mode-pairing at each branch of the mode tree.
%% Minnotte showed that tests on the mode tree away from the bifracation points are
%% less powerful than at the split points.
Minnotte's idea is related to Silverman's \citeyear{silverman81} bootstrap test,
but while Silverman counted the total number of modes in each bootstrap sample,
Minnotte computes and compares the size of the ``excess mass'' in the region
of the two modes that are about to be combined.
This idea works well in one dimension, but
its extension to two or more dimensions will be challenging and
computationally expensive.  A related approach that does not require
bootstrapping was given by \citeN{excessmass91}.  \citeN{marron98}
highlight regions in the mode tree where the derivative of the
kernel estimate is significantly positive or negative to identify
likely modal regions.

We believe the mode tree is a useful framework
in which to attack the clustering problem.  However, the mode tree suffers
from the fact that all possible modes (and hence clusters) are
displayed without a figure of merit.  \citeN{filtered.mode} proposed
a clever enhancement of the mode tree by replacing the kernel estimator
with a locally adaptive MON they call a
filtered kernel estimator.  Their algorithm usually, but not
always, elevates real modes closer to the top of the mode tree.

In the next section, we propose a new
variety of mode tree which is not based on all possible modes
of a kernel estimate, but only those that satisfy a local
optimization problem.


\section{Partial Mixture Fitting and the PDC}

Using probability density functions which are 
not densities is common practice.  
For example, Bayesians employ improper priors,
which integrate to infinity.
Negative probabilities are observed in orthogonal series and
higher-order kernel methods.
Here, the idea is to fit a nonnegative ``density'' which
integrates to a number (usually) less than 1.  Specifically,
consider the partial probability density model,
\begin{equation*}
{\hat f}(x)=w\cdot\phi(\mu,\sigma^2)\,,
\end{equation*}
where $w$ is considered a free parameter on the
interval $(0,1]$.  This model comprises one component of
the full mixture model in Equation (\ref{def:mixnor})
and will be denoted by PDC for {\it partial density component}.

Consider fitting a partial probability
model to the mixture density shown in Figure \ref{mix:pdc:theory}a where $K=3$ and
\begin{equation}
(\mu_k)=\left(-1,0,3\right)^T\qquad
(\sigma_k)=\left(\frac{1}{9},1,\frac{1}{3}\right)^T\qquad
(w_k)=\left(\frac{3}{51},\frac{32}{51},\frac{16}{51}\right)^T.
\label{mix.eg}
\end{equation}
First, observe that no meaningful partial probability model can be
fit by maximum likelihood, because the weight $w$ will
always take the value 1 when the likelihood is
maximized.  This follows since the logarithm of the density
separates the weight $w$ from the other parameters.
Such is not the case, however, if a minimum distance criterion
is employed.  In this paper, integrated squared error (ISE)
or the squared $L_2$-norm is used:
\begin{equation}
  \hbox{ISE}=\int_{-\infty}^{\infty}\left[{\hat f}(x)-f(x)\right]^2dx\,.
  \label{ise:def}
\end{equation}
Integrated squared error is a commonly used
goodness-of-fit measure in nonparametric density estimation.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig03a.eps}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig03b.eps}
   \caption{
	(a)  Mixture of 3 normals.
        (b)  Two partial normal fits according to integrated
	squared error (dashed lines).  Parameter values are
	given in the text.}
    \label{mix:pdc:theory}
 \end{center}
\end{figure}

When $f$ is a normal mixture, a simple closed-form
expression exists for the theoretical ISE, which is of
interest as the solution to a consistent data-based
approach as $n\rightarrow\infty$.
First expand Equation (\ref{ise:def}) to
\begin{equation}
\hbox{ISE}(\mu,\sigma,w)=
	\int\left[ {\hat f}(x)^2 - 2{\hat f}(x)\,f(x) + f(x)^2 \right] dx\,.
\label{ise:expand}
\end{equation}
Recall the identity
\begin{equation*}
\int_{-\infty}^{\infty}\phi(x|\mu_1,\sigma_1^2)\,\phi(x|\mu_2,\sigma_2^2)\,dx =
	\phi(0|\mu_1-\mu_2,\sigma_1^2+\sigma_2^2)\,.
\end{equation*}
Then for a normal PDC, the first integral in Equation (\ref{ise:expand})
equals $w^2\phi(0|0,2\sigma^2)$ or $w^2/(2\,\sigma \sqrt{\pi})$.
The third integral equals
\begin{equation*}
  \int_{-\infty}^{\infty} f(x)^2dx\equiv R(f)=
  \sum_{k=1}^K\sum_{\ell=1}^K w_k\,w_{\ell}\,
  \phi(0|\mu_k-\mu_{\ell},\sigma_k^2+\sigma_{\ell}^2)\,.
\end{equation*}
A similar result exists for the second integral. Combining gives
\begin{equation}
  \hbox{ISE}(\mu,\sigma,w) = \frac{w^2}{2\,\sigma \sqrt{\pi}} -
	2w\sum_{k=1}^K w_k\,\phi(0|\mu-\mu_k,\sigma^2+\sigma_k^2) + R(f)\,.
  \label{ise:theory}
\end{equation}
Since the last term, $R(f)$, does not depend on any of the parameters
of the PDC, $R(f)$ may be ignored when finding the best set of parameters.

With the parameters in Equation (\ref{mix.eg}),
minimization of Equation (\ref{ise:theory}) over $(\mu,\sigma,w)$ 
leads to one of two partial densities:  $(.707,2.066,1.061)$
or $(2.99,.348,.326)$; see Figure \ref{mix:pdc:theory}b.
The second PDC %%% almost exactly
matches the third component values of $(3,.333,.314)$.  The first PDC
is a ``best'' fit to the entire density, but the parameters
differ from the theoretical moments of the mixture
($\mu=.882$ and $\sigma=1.663$).  Also note that for the unconstrained
best $L_2$ fit, $w^*$ turns out to be greater than 1.
There is no local PDC solution
that fits either of the left two bumps alone.
As \citeN{scott.l2e.99} %%% Scott (1999--CS proc)
has noted, PDC's are not attracted to mixture
components with mass below a certain threshold, unless a
component is sufficiently isolated.
Thus fitted PDC's may or may not correspond to individual
components in the underlying mixture density.
%% [Perhaps a bit more....say six sigma.]

Can the PDC be modified to find all three modes in
the mixture?  An affirmative answer lies in recognizing that the
minimization of ISE does not preclude optimizing over
a {\it subset} of the three PDC parameters, $(\mu,\sigma,w)$.
Since greatest interest lies with a cluster's location,
the parameter $\mu$ is always optimized.  By fixing either
$\sigma$ or $w$ or both, can all three components be resolved?
%%% (The question posed is theoretical, that is, $n=\infty$.)
In Figure \ref{pdc:6cases}, partial PDC fits
were computed by minimizing criterion (\ref{ise:theory}),
first with $\sigma$ fixed and optimizing over $(\mu,w)$,
and then with $w$ fixed and optimizing
over $(\mu,\sigma)$, for the same 1,000 random
initial values.  (The sampling ranges for the initial values
of $(\mu_j,\sigma_j,w_j)$, $j=1,\ldots,1000$,
were $(-4,3)$, $(.05,2)$, and $(.05,1)$, respectively.
The sampling distributions were not uniform for
$\sigma$ and $w$, but biased towards smaller values.)

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig04aa.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig04ab.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig04ac.eps} \\
   \includegraphics[scale=0.195,angle=270]{FIGS/fig04ad.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig04ae.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig04af.eps}
   \caption{PDC parameter estimates for 1,000 random
	starting points.  In the top three frames, the
	minimization is with respect to $\mu$ and $w$ only,
	with $\sigma$ fixed.  In the bottom three frames, the
	minimization is with respect to $\mu$ and $\sigma$ only,
	with $w$ fixed.}
   \label{pdc:6cases}
 \end{center}
\end{figure}

In Figure \ref{pdc:6cases}a, which displays
the 1000 points, $(\mu^*_j,\sigma_j)$, all three components are
clearly indicated when $\sigma<0.2$.  When $0.4<\sigma<1$,
the left two components are not separated.  Only
one mean location is found when $\sigma>1.05$.
Observe that the estimates of the mean locations vary smoothly
with $\sigma$, except at several critical values
where the ``branch'' of the curve ends.
The same structure is apparent with the corresponding $(\mu^*_j,w^*_j)$ points shown
in Figure \ref{pdc:6cases}b,
which is easily explained by noting the high correlation between
$(\sigma_j,w^*_j)$ in Figure \ref{pdc:6cases}c for each ``branch.''

A similar but less compelling story is
apparent in Figures \ref{pdc:6cases}(d--f), where the optimization
is over $\mu$ and $\sigma$ with $w$ fixed at its initial random value.
In frame \ref{pdc:6cases}d, only one mean location is found when $\sigma_j^*>0.61$.
The expected branch near 0 appears when $\sigma_j^*<.30$, but
disappears when $\sigma_j^*<.18$.  The final mode at $x=-1$ appears when $\sigma_j^*<0.17$.
The pair $(\sigma^*_j,w_j)$ are still highly correlated
in  \ref{pdc:6cases}f, but the view of $(\mu_j^*,w_j)$ in
\ref{pdc:6cases}e is much less clear as the mode at 0 barely appears at all,
and the ``combined'' mode at $-.26$ persists after the modes at 0 and -1
have appeared.  Empirically, we prefer $(\mu_j^*,\sigma_j,w_j^*)$ to $(\mu_j^*,\sigma_j^*,w_j)$.

Since greatest interest lies with the location parameter, $\mu$,
the simulation was repeated with initial values of
the three parameters chosen at random, but optimization
limited to $\mu$ alone.
In Figure \ref{pdc:2cases}a, the tree-like structure of $(\mu^*_j,\sigma_j)$
is apparent as before.  The 1000 random standard deviation and weight
inputs are shown in Figure \ref{pdc:2cases}c.
A 3-D view of $(\mu^*_j,\sigma_j,w_j)$ confirms
that Figures \ref{pdc:2cases}a and \ref{pdc:6cases}a are the same,
modulo different random starting values.  Apparently,
the initial value of $w$ does not affect the value of $\mu^*$.
A closer look at Equation (\ref{ise:theory}) reveals that
when $\sigma$ is fixed, $w$ and $\mu$ are uncoupled in the second term.
Hence, the best value of $\mu$ is independent of $w$.
In particular, any fixed or random value may be used for
$w$ during the optimization.
Given $\mu^*$, Equation (\ref{ise:theory}) takes the simple form
$c_1w^2-2c_2w$, and $w^*$ is easily seen to equal $c_2/c_1$,
yielding Figure \ref{pdc:6cases}b again.
Numerical optimization over $\mu$ alone
(i.e., without $w$) significantly reduces the computational effort.
%% With practical numerical optimization, the particular one of the several $\mu^*$'s that
%% exist at the given level of $\sigma$ is found by numerical optimization
%% may depend on the initial guess for $w$;  a random start for
%% $w$ (but not part of the nonlinear optimization) or a fixed small
%% choice such as $w=0.2$ may be recommended.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig05a.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig05.5a.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig05b.eps}
   \caption{PDC parameter estimates for 1,000 random
	starting points.  The minimization is with respect
	to $\mu$ only.  (a)  $(\mu^*,\sigma)$.
	(b) $(\mu^*,w)$.  (c) $(\sigma,w)$.}
    \label{pdc:2cases}
 \end{center}
\end{figure}


While each may be useful in practice, after reviewing many figures
such as \ref{pdc:6cases} and \ref{pdc:2cases}, we find that optimizing
over $\mu$ and $w$ with $\sigma$ fixed is preferable.
%% While explicit computation of $w^*$ is not necessary if the desired plot
%% is of $(\mu^*,\sigma)$, we also find the plot of $(\mu^*,w^*)$ is more interpret.
In particular, we find the plot of $(\mu^*,w^*)$ is more interpretable
than $(\mu^*,\sigma)$.
In practice, these diagrams must be
estimated from data, and they will be noisier.
A data-based version of PDC is the next topic.  %%% of the next section.


\section{Stochastic Mode Tree}

The ISE criterion for the PDC
may be estimated with data $\{x_1,\ldots,x_n\}$
from a {\it general} density $f(x)$ by noting that
the second integral in Equation (\ref{ise:expand})
is simply minus two times the expectation $\hbox{E}{\hat f}(X)$
at a random point $X$.
A simple unbiased estimator may be obtained by averaging
${\hat f}(\cdot)$ over the sample data.
As before, the third integral in Equation (\ref{ise:expand})
may be omitted as it is constant with respect to the parameters.
Thus we obtain the following data-based criterion for numerical optimization:
\begin{equation}
  \LTEs(\mu,\sigma,w) = \frac{w^2}{2\,\sigma\sqrt{\pi}} - \frac{2\,w}{n}\sum_{i=1}^n
	\phi(x_i|\mu,\sigma^2)\,.
  \label{sample:l2e}
\end{equation}
The notation $\LTE$ for parametric estimation
was introduced by \citeN{l2e} %% Scott (1998)
to denote the $L_2$ estimation error criterion.
Note that when $\sigma$ is fixed in Equation (\ref{sample:l2e}),
then $w$ and $\mu$ are uncoupled just as in the theoretical
version in Equation (\ref{ise:theory}).  Thus $\mu^*$ may be
computed with $w$ fixed at any value, and then $w^*$ computed
as $c_2/c_1$, using the sample version of $c_2$.

A sample of size $n=1275$ was generated from
the mixture density described in Equation~(\ref{mix.eg}).
Various combinations of the parameters were estimated using random
starts.  %%%  and various subsets of parameters to optimize over.
In Figure \ref{smt:sim.data}, the stochastic mode tree (SMT) for this
simulated dataset is displayed.  While the 3 main modes are clearly evident, so are
a number of small features not in the true density but only
in the sample.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig06a.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig06b.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig06c.eps} \\
   \includegraphics[scale=0.195,angle=270]{FIGS/fig06d.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig06e.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig06f.eps}
    \caption{
	(a--c)  1275 partial density fits $(\mu_j^*,\sigma_j,w_j^*)$ to simulated
		data, with optimization over $\mu$ and $w$ but not $\sigma$:
		(a)  $(\mu_j^*,\sigma_j)$
		(b)  $(\mu_j^*,w_j^*)$
		(c)  The estimates are plotted on the density scale
			at $\left(\mu_j^*,w_j^*\phi(\mu_j^*|\mu_j^*,\sigma_j^2)\right)$.
	(d--f)  1275 partial density fits $(\mu_j^*,\sigma_j^*,w_j)$ with
		the same random starts:
		(d)  $(\mu_j^*,\sigma_j^*)$
		(e)  $(\mu_j^*,w_j)$
		(f)  $\left(\mu_j^*,w_j\phi(\mu_j^*|\mu_j^*,{\sigma_j^*}^2)\right)$.}
    \label{smt:sim.data}
 \end{center}
\end{figure}

Figure \ref{smt.mt:sim.data} displays both the mode tree based on the kernel
estimator as well as the stochastic mode tree.
Note how the stochastic mode tree suppresses
a number of ``potential'' modes in these data.
Specifically, the SMT acts as a partial inferential tool (filter) for gauging
the weight of each mode or cluster of a particular size.  For example, the leftmost
data point at $x=-3.28$ exists as a potential cluster in both figures, but its significance
is obviously nil in the SMT version.  Conversely, the cluster at $x=-1$
is more prominent in the SMT.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig06aa.eps}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig06ab.eps}
    \caption{
	(a)  Stochastic mode tree of simulated data.
	(b)  Original mode tree, but drawn without horizontal connections and
		on a linear scale to facilitate comparison.}
    \label{smt.mt:sim.data}
 \end{center}
\end{figure}



\section{Special Cases}

In this section, we briefly examine the consequences of
using a partial normal component to fit non-Gaussian mixtures.
In the multivariate case, the use of $\Sigma=h\cdot I_d$ is examined
when the local Hessian is not spherical.

A sample of size 1000 was generated from the beta mixture
$\frac{2}{5}\hbox{B}(4,8)+\frac{3}{5}\hbox{B}(10,3)$;  see
Figure \ref{beta:mixture}.
Fitting partial normal components (not the true density)
gives good location estimates for the smaller values of
$\sigma_j$ and $w_j$, i.e., for the more local PDC fits.
%% Again, the optimal weight estimates are
%% less interpretable when the parametric form is not correct.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig08a.eps}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig08b.eps}
    \caption{
	(a)  Beta mixture and sample histogram with $n=1000$.
	(b)  Stochastic mode tree of $(\mu_j^*,w_j^*)$ from
		Normal partial density fits, with arrows indicating the true modal locations.}
    \label{beta:mixture}
 \end{center}
\end{figure}

In the multivariate case, consider a partial density component
model that is Gaussian with covariance matrix proportional to the
identity matrix.  Here it is demonstrated that this model can
find the correct location of Gaussian components for which
the covariance matrix is not proportional to the identity.
For an equal-weight three-component bivariate normal mixture,
a sample of size 1200 was generated and
is shown in Figure \ref{norm2d:mixture}a.
(The mean vectors from left to right are $\binom{-3}{\ \,\,3}$,
$\binom{0}{0}$, and $\binom{\ \,5.0}{-1.5}$.  The correlation
for the first is $.75$ and 0 for the others.  The variances
are all 1 except for $\sigma_y^2=4$ in the last.)
The estimated partial normal components shown in Figure \ref{norm2d:mixture}(b,c)
give good location estimates.  Note that the optimal weight estimates are
less interpretable when the parametric form is not correct.
In particular, fitting $w\cdot\hbox{MN}(\mu,\sigma^2\cdot I_2)$
to data with very large $\sigma$ often results in $w^*>1$.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig09a.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig09b.eps}
   \includegraphics[scale=0.195,angle=270]{FIGS/fig09c.eps}
    \caption{
	(a)  Sample of size 1200 from a bivariate mixture of
		three Normals.
	(b)  Stochastic mode tree of the bivariate Normal
		partial density fits.  A view of
		the points $(\mu_1^*,w^*)$ is shown.
	(c)  A view of the points $(\mu_2^*,w^*)$ is shown.}
    \label{norm2d:mixture}
 \end{center}
\end{figure}

\section{Examples}

\citeN{ripley} uses five measurements on crabs to illustrate
various discrimination algorithms, comparing blue and orange crabs as well as
male and female sexes.  One of the techniques employed is linear
discriminant analysis, which is optimal if the data are
normal, $MN(\mu_k,\Sigma)$.  Here, attention is focused on
the group of male blue crabs ($n=50$). 
The stochastic mode
tree was computed on the blurred and sphered data; see Figure \ref{smt:crabs.bm}.
The figure suggests there may be 2--4 clusters in these data.
The inference of more than one cluster is supported by examining an averaged shifted
histogram (ASH,\citeN{mvde}) of the first two principal components
(explaining 93.8\% of the variance) of the raw data.
In general, the number of clusters will not decrease
as other variables are added. 
(An explanation for the observed ``gap''
in Figure  \ref{smt:crabs.bm}a for $w^*\in(0,.113)$ may be found in
Section \ref{section:infer}.)


\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig10a.eps}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig10b.eps}
    \caption{
	(a)  Projection of the stochastic mode tree of 5-D male blue crab data
		onto the petal width variable.  The horizontal dotted line is at $w^*=0.113$.
	(b)  ASH estimate of the first two principal components of the
		raw male blue crab data.}
    \label{smt:crabs.bm}
 \end{center}
\end{figure}

Figure \ref{smt:iris} displays a portion of the stochastic
mode tree along the petal length variable
for the well-known 4-D Fisher iris data, which were also blurred.
Also shown is a bivariate ASH of that variable and another.  The three
species are indicated, but other modes are apparent.  We return
to this and other examples in the next section.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig11a.eps}
   \includegraphics[scale=0.3,angle=270]{FIGS/fig11b.eps}
    \caption{
	(a)  Portion of the stochastic mode tree of iris data along
		the petal length.
	(b)  ASH estimate of the petal length and width for the
		{\it Setosa}, {\it Versicolor}, and {\it Virginica} species.}
    \label{smt:iris}
 \end{center}
\end{figure}

%%%  mode tree; stochastic mode tree; mixture K=1, K=2, EM and L2E.

%%%   galaxy, stamp, ....


\section{Inferential Possibilities\label{section:infer}}

When $\sigma$ is so small that a data point $x_i$
is ``distant'' from its nearest neighbor, then one PDC solution is
$\mu^*=x_i$;  hence, from equation (\ref{sample:l2e}),
\begin{equation*}
  \LTEs(\mu^*,\sigma,w) = \frac{w^2}{2\,\sigma\sqrt{\pi}} -
	\frac{2\,w}{n\sigma\sqrt{2\pi}}\,.
\end{equation*}
It follows that $w^*=\sqrt{2}/n$, a constant irrespective of
the particular (small) value of $\sigma$.
Thus, if one wished to exclude isolated
clumps of points with fewer than 3 or 4 members, a simple
threshold rule can be developed for excluding any PDC solution
with $w^*\le3\sqrt{2}/n$.  For $\bx\in\Re^d$, a similar argument
shows that
\begin{equation}
w^*=2^{d/2}/n  \label{w.threshold}
\end{equation}
for isolated points.  This formula is quite accurate in practice.
For example, if $d=5$ and $n=50$,
then $w^*=0.113$; cf.~the vertical gap in Figure \ref{smt:crabs.bm}a.
Note that $w^*=0$ is also a stationary point for this problem.
The gap for the Iris data is given by $w^*=2^{4/2}/150=0.027$,
which is (barely) visible in Figure \ref{smt:iris}a.

A more sophisticated line of reasoning explaining
``noise'' in the SMT may be developed.
Suppose a PDC solution, $(\mu^*,\sigma)$,
exists and $\sigma$ is small.
Consider the original $\LTE$ criterion:
\begin{equation}
\int {\hat f}(x)^2dx - 2 \int {\hat f}(x)\,f(x)\,dx\,.
\label{lte.theory.1}
\end{equation}
Taking a Taylor's series about the point $x=\mu^*$
for a Normal PDC, Equation (\ref{lte.theory.1}) becomes
\begin{equation*}
\frac{w^2}{2\sqrt{\pi}\sigma} - 2 w \int_x \phi(x|\mu^*,\sigma^2) \left[
   f(\mu^*) + (x-\mu^*) f'(\mu^*) + \frac{1}{2}(x-\mu^*)^2 f''(\mu^*) + \cdots \right]\,dx
\end{equation*}
or, since $\sigma$ is small
%% , the PDC is symmetric around $x=\mu^*$,
and $\mu^*$ and $\sigma^2$ are the moments of the normal PDC,
\begin{equation*}
\frac{w^2}{2\sqrt{\pi}\sigma} - 2 w \left[
	f(\mu^*) + \frac{1}{2}\sigma^2 f''(\mu^*) + O(\sigma^4) \right]\,.
\end{equation*}
Retaining only the leading terms, the optimal value of the weight approximately satisfies
\begin{equation*}
w^* = 2\sigma\sqrt{\pi}  f(\mu^*) + \sigma^3\sqrt{\pi} f''(\mu^*)\,.
\end{equation*}
Again, this argument may be extended to the multivariate case.
For $\bx\in\Re^d$, let $\hat f$ be a normal PDC with $\Sigma=\sigma^2\,I_d$.
Then a similar argument leads to
\begin{equation*}
w^* = \left(2\,\sigma\sqrt{\pi} \right)^{d} f(\mu^*)+O(\sigma^{d+2})\,.
\end{equation*}


Intuitively, a ``false'' mode in the SMT gives rise
to points $(\mu^*,w^*)$ which approximately satisfy the leading term
$w^* = 2\sigma\sqrt{\pi} f(\mu^*)$.  For true modes, on the other hand,
the linear approximation will not be sufficient.
To examine this hypothesis, we re-examined the SMT results
in Figure \ref{smt.mt:sim.data}.
Using the true density, we computed $f(\mu^*_j)$, then a predicted optimal
weight for each solution by ${\hat w}_j = 2 \sigma \sqrt{\pi} f(\mu^*_j)$.
Using XGobi \cite{xgobi}, we
plotted $(w^*_j,{\hat w}_j)$, and brushed all the points falling
on the 45-degree line near the origin.
In Figure \ref{mix.eg:signif}, the highlighted points as
well as the correspondingly highlighted stochastic mode tree are
displayed.  Clearly, this approach successfully eliminates most of
the background clutter while retaining all and only the 3 components.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.45,angle=0]{FIGS/foo1.eps}
   \includegraphics[scale=0.45,angle=0]{FIGS/foo2.eps}
    \caption{
	(a)  The highlighted (red) points satisfy the local Taylor's Series prediction
		for $w^*$ given $\sigma$, using the true density for $f(\mu^*)$.
	(b)  The corresponding stochastic mode tree for the 1275 points with
		color linked to the left frame.}
    \label{mix.eg:signif}
 \end{center}
\end{figure}

In practice, an estimate of the density is required to
find $\hat w$.  Either a kernel estimator, a mixture estimator,
or a nearest-neighbor estimator may be used for this task.  While
we eventually hope to use a mixture estimate for this task, we
chose a $k$-th nearest-neighbor estimator with fairly small $k$.
While such an estimate is noisier than a kernel estimate, it is
better able to cope with different sized bumps in high dimensions
with only a little tweaking.  

For the (blurred) iris data, we chose $k=8$ and
used XGobi to brush points in the $(w_j^*,{\hat w}_j)$ plot
near the origin on the 45-degree line;  see Figure 
\ref{iris.eg:signif}a.  When the corresponding points in
the stochastic mode tree are plotted in Figure
\ref{iris.eg:signif}b, the three species are clearly
highlighted, together, perhaps, with a hybrid species
of {\it Versicolor} and {\it Virginica} (\citeN{thomp2000}, pp.~250--256).
(Compare to Figure \ref{smt:iris}b.)

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.45,angle=0]{FIGS/iris.bl.infer.eps}
   \includegraphics[scale=0.45,angle=0]{FIGS/iris.bl.smt.eps}
    \caption{
	(a)  For the blurred 4-D Iris data, highlighted (red) points
		approximately satisfy the local
		Taylor's Series prediction for $w^*$ given $\sigma$,
		using an $8$-NN density estimate for $f(\mu^*)$.
	(b)  Corresponding brushed stochastic mode tree projected onto petal length.}
    \label{iris.eg:signif}
 \end{center}
\end{figure}


For the simulated 3-component bivariate example, we computed $\hat f$
using an 8-th nearest-neighbor estimate, rather than the true density.
The brushed $(w^*_j,{\hat w}_j)$ plot and the highlighted stochastic
mode tree are displayed in Figure \ref{sim2d.eg:signif}.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.45,angle=0]{FIGS/sim2d.infer.eps}
   \includegraphics[scale=0.45,angle=0]{FIGS/sim2d.smt.eps}
    \caption{
	(a)  For the 3-component bivariate simulated dataset,
		highlighted (red) points approximately satisfy the local Taylor's Series prediction
		for $w^*$ given $\sigma$, using $8$-NN density estimate for $f(\mu^*)$.
	(b)  The corresponding stochastic mode tree against $\mu_2$.}
    \label{sim2d.eg:signif}
 \end{center}
\end{figure}

Finally, we re-examined the blurred male blue crab data.  We computed $\hat f$
using an 4-th nearest-neighbor estimate.
The brushed $(w^*_j,{\hat w}_j)$ plot and the highlighted stochastic
mode tree are displayed in Figure \ref{crab.eg:signif}.  Note that
all the points at $w^*=0$ fall at the origin in the left frame.
Three groups seem to be supported by this result.

\begin{figure}[htbp]
 \begin{center}
   \includegraphics[scale=0.45,angle=0]{FIGS/crab.infer.eps}
   \includegraphics[scale=0.45,angle=0]{FIGS/crab.smt.eps}
    \caption{
	(a)  Highlighted (red) points satisfy the local Taylor's Series prediction
		for $w^*$ given $\sigma$, using $4$-NN density estimate for $f(\mu^*)$.
	(b)  The corresponding stochastic mode tree against $\mu_2$.}
    \label{crab.eg:signif}
 \end{center}
\end{figure}

Results from other well-known univariate datasets are not shown here,
but are available from our web site.


\section{Conclusions\label{section:conclusions}}

Modes are an excellent summary of data.
A kernel estimator can have many modes, but as the
smoothing parameters change, so do the number and location of the modes.
Determining which of these modes are ``real'' is not an easy task.
The ordinary mode tree displays all possible modal locations without criticism.
In one or two dimensions, some tests and procedures are available.
In more dimensions, mixture models are of particular interest,
but begin to experience practical problems as numerical fitting
produces more singular components.
%%%% gets more and more singular, even with only moderately hard problems.

The stochastic mode tree ``probes'' the (multivariate) data density
surface for normal patches and modes that are located and scaled
to provide best local fits for a given (random) scale.
Alternatively, choosing a polynomial model (e.g.~centered Beta)
for the partial density component may make sense from a Taylor's series' point of view.
Plotting $w^*$ against $\mu^*$ reveals actual modal structure
more clearly.  Questions of whether an individual mode
exists may be tested by fitting mixtures with
components located at the indicated positions (i.e.
don't optimize over the locations);  see \citeN{ss.mix.99}.
Another possible method is to examine mode trees from
bootstrap samples;  see \citeN{mode.forest} for details.
The interactive techniques described in Section \ref{section:infer}
give very promising results as well.

The algorithm as stated is rotationally invariant, not because the data are,
but because the covariance model of the PDC is spherical.
A local estimate of the shape of the density there may be obtained
by computing the negative of the Hessian matrix at each sample mode.
An improved PDC might use that matrix for the covariance,
and iterate several times.  Finally, the estimates could be used to
initialize a MON fit that does not optimize over the mean locations.

The techniques described here are exploratory, as is
cluster analysis generally, so the approach
in Section \ref{section:infer} is quite reasonable.
We also have attempted fitting 2-component PDC's, with random starts.
This line of inquiry is promising but much more work is required to
understand how to extract the meaningful interpretations from
the collection of solutions.

%% We have an open question:  How hard is it to constrain {\it mode} locations
%% versus the mixture {\it mean} locations?

%% One novel use of the stochastic mode tree is feature testing.
%% Consider the possibility that the right mode could actually
%% be two (as suggested in both figures).  Fitting a one-mode and two-mode mixture,
%% and using BIC, we see that there is no support for the
%% extra mode.  Fitting mixtures is much easier if the mixture
%% locations are known.


%% [SIZER stuff of Marron]

%% [Mention PDC's with 2 components.....lots of non-converging solutions,
%% but there are interesting solutions.....rerun with sigma's fixed!!]


%% ***************** re-compute  mt.all   which only has one
%% solution I think**************

%%  For $N(0,1)$ data, compute the height and locations of
%% the top 3 ``branches'' --- using BIC how many are significant?

%%  

\section{Acknowledgments}

This research was supported in part by NSF Grant 99-71797
and was performed while the first author was on sabbatical
at the NSA.

\bibliographystyle{chicago}
\bibliography{stoch.mt}

\section*{Keywords}
%% {\tt Modes, number of clusters, mixture models, minimum distance estimation, kernel density}
{\tt Kernel density; Minimum distance estimation; Mixture models; Modes; Number of clusters}

\end{document}

latex stoch.mt.tex
bibtex stoch.mt
latex stoch.mt.tex
dvips stoch.mt.dvi
ghostview stoch.mt.ps&

latex stoch.mt.tex;bibtex stoch.mt;latex stoch.mt.tex;dvips stoch.mt.dvi
