1

Topological mixture estimation
Steve Huntsman∗

arXiv:1712.04487v1 [stat.ME] 12 Dec 2017

BAE Systems, 4301 North Fairfax Drive, Arlington, Virginia 22203
(Dated: December 14, 2017)
Density functions that represent sample data are often multimodal, i.e. they exhibit more than
one maximum. Typically this behavior is taken to indicate that the underlying data deserves a
more detailed representation as a mixture of densities with individually simpler structure. The
usual specification of a component density is quite restrictive, with log-concave the most general
case considered in the literature, and Gaussian the overwhelmingly typical case. It is also necessary
to determine the number of mixture components a priori, and much art is devoted to this. Here, we
introduce topological mixture estimation, a completely nonparametric and computationally efficient
solution to the one-dimensional problem where mixture components need only be unimodal. We repeatedly perturb the unimodal decomposition of Baryshnikov and Ghrist to produce a topologically
and information-theoretically optimal unimodal mixture. We also detail a smoothing process that
optimally exploits topological persistence of the unimodal category in a natural way when working
directly with sample data. Finally, we illustrate these techniques through examples.

I.

INTRODUCTION

Let D(Rd ) denote a suitable space of continuous probability densities (henceforth merely called densities)
on Rd .
P
d
◦
d M
◦
M
A mixture on R with M components is a pair (π, p) ∈ ∆M × D(R ) , where ∆M := {π ∈ (0, 1] : m πm = 1}; we
write |(π, p)| := M , and note that π cannot have any components equal to zero. The corresponding mixture density
PM
is hπ, pi := m=1 πm pm . The Jensen-Shannon divergence of (π, p) is [15]
J(π, p) := H (hπ, pi) − hπ, H(p)i

(1)

R

where H(p)m := H(pm ) and H(f ) := − f log f dx is the entropy of f .
Now J(π, p) is the mutual information between the random variables Ξ ∼ π and X ∼ hπ, pi. Since mutual
information is always nonnegative, the same is true of J. The concavity of H gives the same result, i.e. H (hπ, pi) ≥
hπ, H(p)i. If M := |(π, p)| > 1 and



πM −1 pM −1 + πM pM
,
(π̂, p̂) := (π1 , . . . , πM −2 , πM −1 + πM ) , p1 , . . . , pM −2 ,
πM −1 + πM
then is easy to show that J(π̂, p̂) ≤ J(π, p).
We say that a density f ∈ D(Rd ) is unimodal if f −1 ([y, ∞)) is either empty or contractible (within itself, or
equivalently, has the homotopy type of a point) for all y: for d = 1, this simply means that any nonempty sets
f −1 ([y, ∞)) are intervals. We call a mixture (π, p) unimodal iff each of the component densities pm is unimodal.
The unimodal category ucat(f ) is the least integer such that there exists a unimodal mixture (π, p) with ucat(f )
components satisfying f = hπ, pi: in this event we write (π, p) |= f . The unimodal category is a topological invariant
that generalizes the geometric category and is related to the Lusternik-Schnirelmann category [4, 12].
The preceding constructions naturally lead us to consider the unimodal Jensen-Shannon divergence
J∩ (f ) :=

sup J(π, p)

(2)

(π,p)|=f

as a simultaneous measure of both the topological and information-theoretical complexity of f , and
(π∩ , p∩ ) := arg max J(π, p)

(3)

(π,p)|=f

as an information-theoretically optimal topological mixture estimate (TME). The natural questions are if such an
estimate exists (is the supremum attained?), is unique, and if so, how to perform TME in practice. In this paper we
address these questions for the case d = 1, and we demonstrate the utility of TME in examples (see Figures 4 and 5).

∗ Electronic

address: steve.huntsman@baesystems.com

December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

2
II.

RELATED WORK

While density estimation enables various clustering techniques [2, 20], mixture estimation is altogether more powerful than clustering: e.g., it is possible to have mixture components that significantly and meaningfully overlap. For
example, a cluster with a bimodal density will usually be considered as arising from two unimodal mixture components that are of interest in their own right. In this light and in view of its totally nonparametric nature, our
approach can be seen as particularly powerful, particularly when coupled with topological density estimation and
deblurring/reblurring (see §VII and §VIII).
Still, even for clustering (even in one dimension, where an optimal solution to k-means can be computed efficiently
[13, 28, 37]), determining the number of clusters in data [26] is as much an art as a science. All of the techniques
we are aware of either require some ad hoc determination to be made, require auxiliary information (e.g., [35]) or are
parametric in at least a limited sense (e.g., [34]). While a parametric approach allows likelihoods and thus various
information criteria [7] or their ilk to be computed for automatically determining the number of clusters, this comes at
the cost of a very strong modeling assumption, and the criteria values themselves are difficult to compare meaningfully
[24].
These shortcomings–including determining the number of mixture components–carry over to the more difficult problem of mixture estimation. [22–24] As an example, an ad hoc and empirically derived unimodal mixture estimation
technique that requires one of a few common functional forms for the mixture components has been recently employed in [25]. Univariate model-based mixtures of skew distributions admit EM-type algorithms and can outperform
Gaussian mixture models [5, 21]: though these generalize to the multivariate case quite effectively (see, e.g., [19]), the
EM-type algorithms are generically vulnerable to becoming trapped in local minima without good initial parameter
values, and they require some model selection criterion to determine the number of mixture components (though the
parameter learning and model selection steps can be integrated as in [10]). A Bayesian nonparametric mixture model
that incorporates many–but not arbitrary–unimodal distributions is considered in [32]. Principled work has been done
on estimating mixtures of log-concave distributions [36]; [8] describes how densities of discrete unimodal mixtures can
be estimated. However, actually estimating generic unimodal mixtures themselves appears to be unaddressed in the
literature, even in one dimension. Indeed, even estimating individual modes and their associated uncertainties or
significances has only been addressed recently [11, 27].

III.

OUTLINE

Given f , the “sweep” algorithm of [4] yields (π, p) |= f . We will show how to repeatedly perturb (π, p) to obtain
(3) using two key lemmas. The first lemma, that J is convex under perturbations of (π, p) that preserve hπ, pi, is
the subject of §IV. The second lemma, a characterization of perturbations of two components of a piecewise affine
and continuous (or piecewise constant) mixture that preserve the predicate (π, p) |= f (i.e., that preserve unimodality
and the mixture density), is the subject of §V. Together, these results entail that greedy unimodality- and densitypreserving local perturbations of pairs of mixture components converge to (3). The proof in §VI is scarcely more
involved than the preceding statement once the lemmas are in hand. Next, in §VII we review the related technique
of topological density estimation (TDE) before showing in §VIII how blurring and deblurring mixture estimates can
usefully couple TDE and TME. Finally, in §IX we produce examples of TME in action.

IV.

CONVEXITY

The following lemma shows that J is convex as we gradually shift part of one mixture component to another.
Lemma. Let |(π, p)| = 3. The function gπ,p : [0, 1] → [0, ∞) defined by



π1 p1 + (1 − t)π2 p2 tπ2 p2 + π3 p3
gπ,p (t) := J (π1 + (1 − t)π2 , tπ2 + π3 ) ,
,
π1 + (1 − t)π2
tπ2 + π3

(4)

satisfies
gπ,p (t) ≤ t · gπ,p (1) + (1 − t) · gπ,p (0).

December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

(5)

3
Proof. The statement as well as the proof of the lemma can be simplified with a bit of extra notation. Define
π12,t := π1 + (1 − t)π2 ;
π23,t := tπ2 + π3 ;
π1 p1 + (1 − t)π2 p2
p12,t :=
;
π12,t
tπ2 p2 + π3 p3
.
p23,t :=
π23,t
Now h(π12,t , π23,t ) , (p12,t , p23,t )i = hπ, pi, so
gπ,p (t) = J ((π12,t , π23,t ) , (p12,t , p23,t ))
= H(h(π12,t , π23,t ) , (p12,t , p23,t )i) − h(π12,t , π23,t ) , (H(p12,t ), H(p23,t ))i
= H(hπ, pi) − h(π12,t , π23,t ) , (H(p12,t ), H(p23,t ))i .
Furthermore, if we write π12 := π1 + π2 , π23 := π2 + π3 , p12 :=

π1 p1 +π2 p2
,
π12

and p23 :=

π2 p2 +π3 p3
,
π23

(6)

then

π12,t = tπ1 + (1 − t)π12 ;
π23,t = tπ23 + (1 − t)π3 ;
tπ1 p1 + (1 − t)π12 p12
p12,t =
;
π12,t
tπ23 p23 + (1 − t)π3 p3
.
p23,t =
π23,t
It is well known that H is a concave functional: from this it follows that
tπ1
(1 − t)π12
H(p1 ) +
H(p12 );
π12,t
π12,t
tπ23
(1 − t)π3
H(p23,t ) ≥
H(p23 ) +
H(p3 ).
π23,t
π23,t
H(p12,t ) ≥

Therefore
gπ,p (t) =
≤
=
=

H(hπ, pi) − h(π12,t , π23,t ) , (H(p12,t ), H(p23,t ))i
H(hπ, pi) − tπ1 H(p1 ) − (1 − t)π12 H(p12 ) − tπ23 H(p23 ) − (1 − t)π3 H(p3 )
t [H(hπ, pi) − h(π1 , π23 ), (H(p1 ), H(p23 ))i] + (1 − t) [H(hπ, pi) − h(π12 , π3 ), (H(p12 ), H(p3 ))i]
t · gπ,p (1) + (1 − t) · gπ,p (0)

as claimed. 
V.

PRESERVING UNIMODALITY

Suppose that (π, p) is a unimodal mixture on R with |(π, p)| > 1. We would like to determine how we can perturb
two components of this mixture so that the result is still unimodal and yields the same density. In the event that
the mixture is piecewise affine and continuous (or piecewise constant) the space of permissible perturbations can be
characterized by the following
Lemma. For 0 ≤ k ≤ N , let yk ∈ [0, ∞) be such that there are integers `, u satisfying 0 < ` ≤ u < N and
0 = y0 ≤ · · · ≤ y`−1 < y` = · · · = yu > yu+1 ≥ · · · ≥ yN = 0.

(7)

(That is, y1 , . . . , yN −1 is a nonnegative, nontrivial unimodal sequence.) Then for 1 ≤ r ≤ N − 1 and ε−
r ≥ 0,
−
yk − δkr ε−
r is nonnegative and unimodal ⇔ εr ≤ yr − (yr−1 ∧ yr+1 ).

December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

(8)

4
Similarly, for ε+
r ≥ 0,
(
+
yk + δkr ε+
r is nonnegative and unimodal ⇔ εr ≤

∞
(yr−1 ∨ yr+1 ) − yr

if ` − 1 ≤ k ≤ u + 1
otherwise.

(9)

Proof (sketch). We first sketch ((8),⇐). Nonegativity follows from 0 ≤ yr−1 ∧ yr+1 ≤ yr − ε−
r . Unimodality
follows from a series of trivial checks for the cases 1 ≤ r < ` − 1, r = ` − 1, r = `, and ` < r < u: the remaining cases
r = u, r = u + 1, and u + 1 < r ≤ N − 1 follow from symmetry. For example, in the case 1 ≤ r < ` − 1, we only need
to show that yr−1 ≤ yr − ε−
r ≤ yr+1 .
A sketch of ((8),⇒) amounts to using the same cases and symmetry argument to perform equally trivial checks.
For example, in the case 1 ≤ r < ` − 1, we have ε−
r ≤ yr − yr−1 ≤ yr − (yr−1 ∧ yr+1 ).
The proof of (9) is mostly similar to that of (8): the key difference here is that any argument adjacent to or at
a point where the maximum is attained can have its value increased arbitrarily without affecting unimodality (or
nonnegativity). 
The example in figure 1 is probably more illuminating than filling in the details of the proof sketch above.

FIG. 1: A unimodal sequence y with unimodality-saturating perturbations ε± respectively indicated above and below.

VI.

ALGORITHM

Theorem. Let −∞ = x−1 < x0 < · · · < xN < xN +1 = ∞ and f be piecewise constant (or affine) over each
[xk , xk+1 ]. Then there is an efficient algorithm to compute (3).
Proof. First, compute (π, p) |= f via the “sweep” algorithm of [4]. Then, for each k ∈ [N ] and pair of mixture
components, compute the largest possible local perturbation of the mixture according to the lemma in §V. By the
lemma in §IV, choosing a resulting maximal value of J and repeating the perturbation-evaluation procedure gives the
desired result in O(M N ) iterations. This result is unique by convexity. 
VII.

TOPOLOGICAL DENSITY ESTIMATION

The obvious situation of practical interest for TME is that a density has been obtained from a preliminary estimation
process involving some sample data. There is a natural approach to this preliminary estimation process called
topological density estimation (TDE) [16] that naturally dovetails with (and as it happens, inspired) TME.
We recall the basic idea here. Given a kernel K and sample data Xj for 1 ≤ j ≤ n, and for each proposed bandwidth
h, we compute the kernel density estimate [9, 33]
n

1X
fˆh;X (x) :=
KXj ,h
n j=1

December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

(10)

5
where Kµ,σ (x) :=

x−µ
1
σ K( σ ).

Next, we compute
uX (h) := ucat(fˆh;X )

(11)

and estimate the unimodal category of the PDF that X is sampled from via
m̂X := arg max µ(u−1
X (m))
m

(12)

where µ denotes an appropriate measure (nominally counting measure or the pushforward of Lebesgue measure under
the transformation h 7→ 1/h).
That is, (12) gives the most prevalent (and usually in practice, also topologically persistent [12, 31]) value of the
unimodal category, i.e., this is a topologically robust estimate of the number of components required to produce the
PDF that X is sampled from as a mixture. While any element of u−1
X (m̂X ) is a bandwidth consistent with the
estimate (12), we typically make the more detailed nominal specification
ĥX := medianµ (u−1
X (m̂X )).

(13)

TDE turns out to be very computationally efficient relative to the traditional technique of cross-validation (CV).
On highly multimodal densities, TDE is competitive or at least reasonably performing relative to CV and other
nonparametric density estimation approaches with respect to traditional statistical evaluation criteria. Moreover,
TDE outperforms other approaches when qualitative criteria such as the number of local maxima and the unimodal
category itself are considered (see Figures 2 and 3). In practice, such qualitative criteria are generally of paramount
importance. For example, precisely estimating the shape of a density is generally less important than determining if
it has two or more clearly separable modes.
j
As an illustration, consider µ(j, m) := m+1
, σ(k, m) := 2−(k+2) (m + 1)−2 and the family of distributions
m

fkm :=

1 X
Kµ(j,m),σ(k,m)
m j=1

(14)

for 1 ≤ k ≤ 3 and 1 ≤ m ≤ 10, and where here K is the standard Gaussian density: see Figure 2. Exhaustive details
relating to the evaluation of TDE on this family and other densities are presented in the software package and test
suite [3]: here, we content ourselves with the performance data for (14) shown in Figure 3.

FIG. 2: The densities fkm in (14) for 1 ≤ k ≤ 3 and 1 ≤ m ≤ 6 over the interval [−0.5, 1.5]. Rows are indexed by k and
columns are indexed by m; the upper left panel shows f11 and the lower right panel shows f36 .

TDE has the very useful feature (shared by essentially no high-performing density estimation technique other than
CV) that it requires no free parameters or assumptions. Indeed, TDE can be used to evaluate its own suitability: for
unimodal distributions, it is clearly not an ideal choice–but it is good at detecting this situation in the first place.
Furthermore, TDE is very efficient computationally.
In situtations of practical interest, it is tempting to couple TDE and TME in the obvious way: i.e., perform them
sequentially and indepdently. This yields a completely nonparametric estimate of a mixture from sample data alone.
However, there is a better way to couple these techniques, as we shall see in the sequel.
VIII.

BLURRING AND DEBLURRING
A.

Blurring

Recall that a log-concave function is unimodal, and moreover that a function is log-concave iff its convolutions with
unimodal functions are identically unimodal. [6, 17, 18] This observation naturally leads to the following question:
December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

6

FIG. 3: (L) Pseudotransparency plots of performance measures relating to the unimodal category for the family (14) depicted
in Figure 2 for n = 500 sample points using a Gaussian kernel. Data for CV is shown in red, while data for TDE is shown in
blue. Each panel represents data on a particular performance measure (left, ucat for TDE alone; middle, ucat for CV alone;
right, the empirical probability that the estimate of ucat is correct) and value of k (top, k = 1; middle, k = 2; bottom, k = 3),
with m varying from 1 to 10. (R) As in the left panel, but for the number of local maxima instead of the unimodal category.

if (π, p) |= f , how good of an approximation to the δ distribution must a log-concave density g be in order to have
(π, p ∗ g) |= f ∗ g? In particular, suppose that g is a Gaussian density: what bandwidth must it have? An answer to
this question of how much blurring a minimal unimodal mixture model can sustain defines a topological scale (viz.,
the persistence of the unimodal category under blurring) that we proceed to illustrate in light of TDE.
In this paragraph we assume that K is the standard Gaussian density, so that Kµ,h ∗ Kµ0 ,h0 = Kµ+µ0 ,(h2 +h02 )1/2 and
fˆh;X ∗ K0,h0 = fˆ(h2 +h02 )1/2 ;X . Write ĥX for the bandwidth obtained via TDE, whether via the nominal specification
−1
ˆ
(13) or any other: by construction we have that inf u−1
X (m̂X ) ≤ ĥX ≤ sup uX (m̂X ). Now if (π, p) |= fĥX ;X , then
m̂X = uX (ĥX ) = |(π, p)|. In order to have (π, p ∗ K0,h0 ) |= f ∗ K0,h0 , it must be that m̂X = uX (ĥX ) = |(π, p)| =
|(π, p ∗ K0,h0 )| = uX ((ĥ2X + h02 )1/2 ), i.e.,
h0 ≤

1/2

2
sup u−1
− ĥ2X
.
X (m̂X )

In particular, we have the weaker inequality involving a purely topological scale:

2 
2 1/2
−1
h0 ≤ sup u−1
(
m̂
)
−
inf
u
(
m̂
)
.
X
X
X
X

(15)

(16)

The preceding considerations generalize straightforwardly if we define uf (h) := ucat(f ∗ K0,h ), where once again K
is a generic kernel. This generalizes (11) so long as we associate sample data with a uniform average of δ distributions.
Under reasonable conditions, we can write uf (0) = ucat(f ), and it is easy to see that the analogous bound is
h0 ≤ sup u−1
f (uf (0)).

(17)

Of course, (17) merely restates the triviality that the blurred mixture ceases to be minimal precisely when the
number of mixture components exceeds the unimodal category of the mixture density. Meanwhile, the special case
furnished by TDE with the standard Gaussian kernel affords sufficient structure for a slightly less trivial statement.
B.

Deblurring/reblurring

The considerations of §VIII A suggest how to couple TDE and TME in a much more effective way than performing
them sequentially and independently. The idea is to use a Gaussian kernel and instead of (13), pick the bandwidth
(−)

ĥX := infµ (u−1
X (m̂X ))
December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

(18)

7
and then perform TME; finally, convolve the results with K0,∆h where

∆h :=

ĥ2X

h
i2 1/2
(−)
− ĥX
.

(19)

This preserves the result of TDE while giving a smoother, less artificial, and more practically useful mixture estimate
than the information theoretically optimal result.
Of course, a similar tactic can be performed directly on a density f by considering its Fourier deconvolution
F −1 (Ff /FK0,h0 ), where F denotes the Fourier transform and h0 is as in (17): however, any a priori justification for
such a tactic is necessarily context-dependent in general, and our experience suggests that its implementation would
be delicate and/or prone to aliasing. Nevertheless, this would be particularly desirable in the context of heavy-tailed
distributions, where kernel density estimation requires much larger sample sizes in order to achieve acceptable results.
In this context it would also be worth considering the use of a symmetric stable density [29, 38] (e.g., a Cauchy
density) as a kernel with the aim of recapturing the essence of (19).

IX.

EXAMPLES

We present two phenomenologically illustrative examples. First, in Figure 4 we consider the n = 272 waiting times
between eruptions of the Old Faithful geyser from the data set in [14]. Then, in Figure 5 we consider the n = 2107
Sloan Digitial Sky Survey g − r color indices accessed from the VizieR database [30] at http://cdsarc.u-strasbg.
fr/viz-bin/Cat?J/ApJ/700/523 and discussed in [1].

FIG. 4: TME applied to n = 272 waiting times between eruptions of the Old Faithful geyser. Left panels, from top to bottom:
area plots of unimodal decompositions obtained by (top) the “sweep algorithm” of [4] on a topological density estimate with
bandwidth given by (13); (second from top) the result of (3) on a topological density estimate with bandwidth given by (13);
(second from bottom) the result of (3) on a topological density estimate with bandwidth given by (18); (bottom) the convolution
of the mixture in the preceding panel with a Gaussian kernel with bandwidth given by (19). Note that three of the four mixture
estimates have the same density, but that the deblurred density is different (this is why we have used the “*” annotation for
this case). Right panels: line plots of the same decompositions.

As suggested, several phenomena are readily apparent from these examples. First, mixtures obtained via the sweep
algorithm are manifestly parity-dependent, i.e., the direction of sweeping matters; second, mixtures obtained via
TME alone exhibit artificial anti-overlapping behavior; third, deblurring followed by reblurring preserves unimodality,
the overall density, a topologically persistent invariant (viz., the unimodal category) and the spirit of informationtheoretical optimality while producing an obviously better behaved mixture; fourth and finally, the various techniques
involved here can significantly shift classification/decision boundaries based on the dominance of various mixture
components.
While the data in Figure 4 is at least qualitatively approximated by a two-component Gaussian mixture, it is clear
that a three-component Gaussian mixture cannot capture the highly oscillatory behavior of the density in Figure 5.
Indeed, this example illustrates how such oscillatory behavior can actually arise from a unimodal mixture with many
fewer components than might naively appear to be required.

December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

8

FIG. 5: TME applied to n = 2107 g − r color indices from [1]. Panels are otherwise as in Figure 4.

Figure 6 shows that there are strong and independent grounds to conclude that the color index data of Figure
5 is produced by a unimodal mixture of three components, with componentwise modes as suggested by TME, and
furthermore that CV undersmooths this data. For each density estimate shown, each componentwise extremum of
the corresponding mixture (3) is virtually identical to one of the local extrema of the density estimate itself: this is a
consequence of the anti-overlapping tendency described above.
In particular, the fourth panel of Figure 6 illustrates that it is possible and potentially advantageous to use TME
as an alternative mode identification technique in the LPMode algorithm of [27]. Furthermore, while we have not
implemented a reliable Fourier deconvolution/reblurring algorithm of the sort hinted at in §VIII B, the fifth and sixth
panels of Figure 6 suggest that this is not particularly important for the narrow task of mode finding/bump hunting.

Acknowledgments

The author thanks Adelchi Azzalini, Robert Ghrist, Subhadeep Mukhopadhyay, and John Nolan for their helpful and
patient discussions. This material is based upon work supported by the Defense Advanced Research Projects Agency
(DARPA) and the Air Force Research Laboratory (AFRL). Any opinions, findings and conclusions or recommendations
expressed in this material are those of the author(s) and do not necessarily reflect the views of DARPA or AFRL.

[1] An, D. et al. “Galactic globular and open clusters in the Sloan Digital Sky Survey. II. Test of theoretical stellar isochrones.”
Astrophys. J. 700, 523 (2009).
[2] Azzalini, A. and Menardi, G. “Clustering via nonparametric density estimation: the R package pdfCluster”. J. Stat. Soft.
57 (2014).
[3] BAE Systems. MATLAB implementation of TDE, available at https://git.io/vyC4l (2017).
[4] Baryshnikov, Y. and Ghrist, R. “Unimodal category and topological statistics.” NOLTA (2011).
[5] Basso, R. M. et al. “Robust mixture modeling based on scale mixtures of skew-normal distributions.” Comp. Stat. Data
Anal. 54, 2926 (2010).
[6] Bertin, E. M. J., Cuculescu, I., and Theodorescu, R. Unimodality of Probability Measures. Kluwer (1997).
[7] Burnham, K. P. and Anderson, D. R. Model Selection and Multimodal Inference: A Practical Information-Theoretic
Approach. 2nd ed. Springer (2003).
[8] Chan, S.-O., et al. “Learning mixtures of structured distributions over discrete domains.” SODA (2013).
[9] Chen, Y.-C. “A tutorial on kernel density estimation and recent advances.” arXiv:1704.03924 (2017).
[10] Figueiredo, M. A. T. and Jain, A. K. “Unsupervised learning of finite mixture models.” IEEE Trans. Pattern Anal. Mach.
Intell. 24, 381 (2002).
[11] Genovese, C. R. et al. “Non-parametric inference for density modes.” J. Royal Stat. Soc. B 78, 99 (2016).
[12] Ghrist, R. Elementary Applied Topology. CreateSpace (2014).
[13] Grønlund, A., et al. “Fast exact k-means, k-medians, and Bregman divergence clustering in 1D.” arXiv:1701.07204 (2017).
[14] Härdle, W. Smoothing Techniques with Implementation in S. Springer (1991).

December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

9

FIG. 6: Structural analyses of the color index data from Figure 5. Top panel: density estimates obtained via CV, the default
LPMode algorithm of [27] (kindly provided by its author), and TDE. Second panel from top: the local extrema of the density
estimates above. Third through fifth panels from top: the result of (3) on the CV, LPMode, and topological density estimates,
respectively, with extrema of components indicated. Bottom panel: the reblurred mixture from the bottom left panel of Figure
5, with extrema of components additionally indicated.

[15] Briët, J. and Harremoës, P. “Properties of classical and quantum Jensen-Shannon divergence.” Phys. Rev. A 79, 052311
(2009).
[16] Huntsman, S. “Topological density estimation.” NOLTA (2017).
[17] Ibragimov, I. A. “On the composition of unimodal distributions.” Th. Prob. Appl. 1, 255 (1956).
[18] Keilson, J. and Gerber, H. “Some results for discrete unimodality.” J. Am. Stat. Ass. 66, 386 (1971).
[19] Lee, S. X. and McLachlan, G. J. “Finite mixtures of canonical fundamental skew-t distributions.” Stat. Comput. 26, 573
(2016).
[20] Li, J., Ray, S., and Lindsay, B. G. “A nonparametric statistical approach to clustering via mode identification.” J. Machine
Learning Res. 8, 1687 (2007).
[21] Lin, T. I., Lee, J. C., and Hsieh, W. J. “Robust mixture modeling using the skew t distribution.” Stat. Comput. 17, 81
(2007).
[22] McLachlan, G. and Peel, D. Finite Mixture Models. Wiley (2000).
[23] McLachlan, G. J. and Rathnayake, S. “On the number of components in a Gaussian mixture model.” WIREs Data Mining
Knowl. Disc. 4, 341 (2014).
[24] Melnykov, V. and Maitra, R. “Finite mixture models and model-based clustering.” Stat. Surv. 4, 80 (2010).
[25] Mints, A. and Hekker, S. “A Unified tool to estimate Distances, Ages, and Masses (UniDAM) from spectrophotometric
data.” Astro. Astrophys. 604, A108 (2017).
[26] Mirkin, B. “Choosing the number of clusters.” Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery 1,
252 (2011).
[27] Mukhopadhyay, S. “Large-scale mode identification and data-driven sciences.” Electronic J. Stat. 11, 215 (2017).
[28] Nielsen, F. and Nock, R. “Optimal interval clustering: application to Bregman clustering and statistical mixture learning.”
IEEE Sig. Proc. Lett. 21, 1289 (2014).
[29] Nolan, J. P. Stable Distributions: Models for Heavy Tailed Data. Birkhaüser (2018).

December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

10
[30] Ochsenbein, F., Bauer, P. and Marcout, J. “The VizieR database of astronomical catalogues.” Astro. Astrophys. Supp.
Ser. 143, 23 (2000).
[31] Oudot, S. Y. Persistence Theory: From Quiver Representations to Data Analysis. AMS (2015).
[32] Rodrı́guez, C. E. and Walker, S. G. “Univariate Bayesian nonparametric mixture modeling with unimodal kernels.” Stat.
Comp. 24, 35 (2014).
[33] Silverman, B.W. Density Estimation for Statistics and Data Analysis. Chapman & Hall (1986).
[34] Sugar, C. A. and James, G. M. “Finding the number of clusters in a dataset: an information-theoretic approach.” J. Am.
Stat. Ass. 98, 750 (2003).
[35] Tibshirani, R., Walther, G., and Hastie, T. “Estimating the number of clusters in a data set via the gap statistic.” J. Roy.
Stat. Soc. B 63, 411 (2001).
[36] Walther, G. “Inference and modeling with log-concave distributions.” Stat. Sci. 24, 319 (2009).
[37] Wang, H. and Song, M. “Ckmeans.1d.dp: optimal k-means clustering in one dimension by dynamic programming.” R J.
3, 29 (2011).
[38] Zolotarev, V. M. and Uchaikin, V. V. Chance and Stability: Stable Distributions and Their Applications. De Gruyter
(1999).

December 14, 2017

Approved for public release; unlimited distribution.
c 2017 BAE Systems. All rights reserved.

