Selective Clustering Annotated using Modes of Projections
Evan Greene Greg Finak Raphael Gottardo


Revisiting the effect of spatial resolution on information content based on classification results
    
M. G. PalacioaCONTACT M. G.Palacio. Email: gpalacio@exa.unrc.edu.ar, S. B. Ferreroa and A. C. Freryb
aUniversidad Nacional de Río Cuarto, Córdoba, Argentina; 
 bUniversidade Federal de Alagoas, Maceió, AL,Brazil

    December 30, 2023
=========================================================================================================================================================================================================================







§ INTRODUCTION


This paper introduces a new algorithm to cluster data in ℝ^p called selective clustering annotated using modes of projections (SCAMP).
Clustering is a common task in data analysis.
Numerous approaches to clustering have been developed over the past century: see 
<cit.>, <cit.>, and <cit.> for 
descriptions of many such methods.
We begin by giving an overview of the SCAMP algorithm, in which we relate it to prior 
work in the clustering literature.

The SCAMP algorithm is inspired by both the mixture modeling 
and density based approach to clustering.
Suppose the rows x_i = (x_i,1,…,x_i,p) of a data matrix X_k× p  
are samples from a distribution F with distribution function (df)

    F(x_i) = ∑_m=1^g π_m F_m(x_i)  with  0 ≤π_m ≤ 1  ∀  m , and ∑_m=1^g π_m = 1 .

In mixture modeling, each component df F_m of the mixture is typically associated 
with an observed cluster.
To cluster the data, a parametric model, such as the multivariate normal 
[<cit.>], uniform normal [<cit.>], multivariate t 
[<cit.>],
or a skew variant of one of these distributions [<cit.>], 
is typically assumed to describe the component distributions F_m.
Non-parametric approaches are also possible: <cit.> suppose the 
mixture is made up of non-parametric components and estimates them with Gaussian 
kernels;
<cit.> and <cit.> take a Bayesian approach to 
estimating mixture models based on unimodal distributions; 
<cit.> take a copula-based approach that influenced the model formulation 
used in this work.
The EM algorithm of <cit.> (or one of its extensions) is then used 
to find local maximizers of the likelihood for a fixed number of clusters g.
The number of components g^* can be determined in terms of the model by picking the g^* that optimizes some criterion (such as BIC) over a user-specified range of g, and subsequently refined by merging clusters with, for example, the method of <cit.> or <cit.>.
The review by <cit.> provides a detailed overview of this approach and references to many recent developments.

SCAMP takes a perspective similar to <cit.>: clusters are defined in terms of unimodality.
In many scientific contexts, unimodality of observations along a measurement coordinate of the data matrix X reflects physical homogeneity of interest.

Using Sklar's theorem, (see Theorem 2.10.9, <cit.>, as well as  <cit.>), 
we can formalize a data model: we assume an observation x is sampled from the mixture with df

    F(x_i) = ∑_m=1^g π_m C_m(F_m,1(x_i,1),F_m,2(x_i,2),…,F_m,p(x_i,p)) ,

where each C_m is a p-copula for the m^th component  subject to the constraint that F_m,j(x_i,j) 
is unimodal for all 1 ≤ m ≤ g and 1 ≤ j ≤ p.

We will reference this formalization (<ref>) as we discuss SCAMP's design choices.
However, SCAMP does not attempt to estimate model parameters using an EM-based approach.
Indeed, without further assumptions on the copulas C_m, using the mixture F to cluster the data is a non-parametric problem.
A standard approach to non-parametric clustering is to estimate the level-sets of F at some level λ > 0 of the mixture df F(x) (see figure <ref>).
Clustering is then achieved by associating observations with local modes.

<cit.> provides an early discussion of how level sets can be used in clustering.
In it, Hartigan noted that different choices of λ induce a hierarchical structure, now commonly called a cluster tree.
He later investigated the problem of  estimating level sets in the two-dimensional case, with an aim of testing for bimodality [<cit.>].
<cit.> also studied this testing problem, leading to the development of the excess mass test.
<cit.> studied theoretical properties of the excess mass, and later the silhouette.
<cit.> developed a method called runt pruning, which aims to construct a pruned cluster tree using a nearest neighbor density estimate.
<cit.> later generalized the pruning method in a graph based clustering approach.
Analysis of the stability of both level set of cluster tree estimates are discussed in [<cit.>].
Two consistent procedures, the first which generalizes single-linkage, the second based on the k-nearest neighbor graph, are analyzed in [<cit.>].
More recently in <cit.>, methods for constructing confidence sets for level sets by <cit.> and <cit.> 
are compared to a new bootstrap based approach.

The level-set approach to clustering makes no distributional assumptions and generalizes naturally to arbitrary dimension.
Compared to mixture modeling, this provides robustness in the event the component distributions of the mixture are specified incorrectly.
On the other hand, approaches that rely on non-parametric density estimators suffer the curse of dimensionality [<cit.>].
For example, <cit.> observe that it is infeasible to use a plug-in binned density-estimator to estimate a cluster tree for a data matrix with ten features since “ten bins per variable in ten dimension would result in 10^10 bins”.
Bandwidth selection also affects the modal structure of the estimated density, which adds an additional complication to non-parametric modal clustering.

The SCAMP algorithm attempts to cluster a data matrix by associating observations with 
components of the mixture distribution (<ref>).   
In deriving these associations, it tries to reap some of the robustness benefits of 
non-parametric modal clustering without suffering excessively from the curse of 
dimensionality. Approaches that combine aspects of non-parametric modal 
clustering and mixture modeling have been proposed by <cit.>,
<cit.>, <cit.>, and <cit.>.
Broadly speaking, these methods rely on fitting a mixture model (typically Gaussian)
for a fixed number of  components, and then using density considerations 
to merge components.

Here, SCAMP takes a partitional approach 
that rests on a severe simplification: the SCAMP algorithm only uses     
one-dimensional density estimates when clustering a dataset.
This simplification carries large costs of its own. 
Most notably, higher-dimensional modal structures that are not axis-aligned cannot 
be detected by SCAMP directly.
However, we argue that the benefits of this simplification outweigh the costs, since it
gives rise to many practical consequences.
For example, SCAMP can cluster datasets with a large number of observations without
sub-sampling, since O(n) implementations exist for many of its constituent
algorithms. We will also argue that SCAMP's tuning parameters are relatively easy 
to understand, as they describe characteristics of one-dimensional densities.
In addition, we note SCAMP can be combined with other methods, such as PCA, to find
structures that are not initially evident along the coordinates of
the data matrix. This is discussed in section <ref> in more depth. 

Underlying the SCAMP algorithm is the assumption that observed data have been generated
by a specific realization of the mixture distribution (<ref>).
Its clustering strategy is based on the following observation:  in a large data matrix X, samples from a component of (<ref>)
with a large mixing weight π_m will produce a sub-collection of rows of X with unimodal empirical distributions along each coordinate projection.
So, to cluster the data matrix, SCAMP begins by searching for such sub-collections of rows of X with unimodal empirical distributions along each coordinate.
We will call the unimodal sub-collections sought by SCAMP α-m-clusters. We pause to define them:


Denote a random subset of the indices 1,2,…,n by 

    ℐ_m^n ≡{(i_(1),…,i_(j))    m ≤ j ≤ n  and  1 ≤ i_(1) < i_(2) < … < i_(j)≤  n} .

When both m and n are understood, we write ℐ in place of ℐ_m^n.



  An α-m-cluster of the matrix X_n× p is any sub-matrix X_ℐ_m^n× p≡ X_ℐ× p such that each of the p coordinate projections of the sub-matrix
  have p-values greater than α when the dip test statistic of <cit.> is computed.


This definition assumes familiarity with the dip test.
A detailed discussion of the dip test occurs later in section <ref>, 
and so we defer elaboration until then.
We note here only that the methods of <cit.> and 
<cit.> both rely on the dip test in their cluster merging methods.

The search for α-m-clusters starts in a similar fashion to that taken by <cit.>:
each coordinate of X is tested for multimodality using the dip test.
If the dip test p-value for a coordinate is smaller than the user selected α, SCAMP then splits that coordinate into modal sub-collections.
In most cases, SCAMP uses the taut string [<cit.>] to determine the split points (section <ref> provides details on how split points are determined in small samples).
By using the dip to guide the search for α-m-clusters, SCAMP's decision to split a coordinate into sub-collections has no dependence on bandwidth: 
so as long as the null hypothesis of unimodality is rejected at the user specified level α, the search for α-m-clusters will continue.

The procedure continues recursively on each modal sub-collection of each multimodal coordinate, until SCAMP encounters a sub-collection of observations that meet one of two stopping conditions:
either the sub-collection satisfies the definition of an α-m-cluster, or the sub-collection contains fewer than m observations.
All branches of a tree that terminate at a leaf with fewer than m observations are pruned.
Notice that after the first recursive step, the density estimates produced by SCAMP are conditional densities that depend on the sequence of modal groupings leading from the root of the data matrix to the specific sub-collection of observations under consideration.
Notice too that a single coordinate of the data matrix can appear multiple times in this search, so long as the multimodality is detected at a particular depth.
In practice we usually constrain the search to split a coordinate only once along any path from a leaf to the root of a tree, both to facilitate cluster interpretation and reduce computational cost.
 
This search generates a forest we call the annotation forest. Each tree in the annotation forest satisfies the following conditions:

  
  * Each tree starts by splitting a coordinate of the data matrix that exhibits multimodality at level α in the root population.
  
  * Each leaf of each tree in the annotation forest is an α-m-cluster.
  
  * Each tree in the forest can be used to define partitional clustering of the data matrix by viewing the leaves of a given tree as the clusters, 
      and combining any observations corresponding to pruned leaves together into a single residual cluster.
  
  * Each observation in the data matrix appears in at most one leaf of each tree of the annotation forest.

Condition <ref> can be interpreted as stating that the SCAMP search defines multiple hierarchical clusterings of the data matrix.
Because of this, we call the trees that make up the annotation forest partition trees.
Unlike agglomerative methods such as single-linkage clustering 
[<cit.>], SCAMP neither requires the specification of nor 
the computation of pair-wise similarities.
While SCAMP has some traits in common with density-based hierarchical clustering approaches such as
DBSCAN [<cit.>], OptiGrid [<cit.>], and 
OPTICS [<cit.>],
we note that the SCAMP search does not require the selection of a metric (nor a distance threshold
ϵ),
does not attempt to order its search according to a best split,
and its stopping condition has no dependence on a bandwidth parameter. 
We refer the reader to the discussion on DBSCAN and OPTICS given by <cit.>.

We will call the α-m-clusters appearing as leaves in this forest as candidate clusters.
Recalling the underlying data model (<ref>),  we view each candidate cluster as a sample from some component C_m of the mixture df F.
However, condition <ref> implies SCAMP must select a subset of the candidate clusters to partition the dataset.
This is because, by construction, a single row of the data matrix appears in multiple candidate clusters when the annotation forest contains more than one partition tree.
According to our data model (<ref>), 
an observation belonging to multiple candidate clusters would indicate it is a sample from multiple different components of the mixture.
Since we do not allow this, SCAMP is forced to make a selection.

To make this selection, we define a preference function and use it to score each candidate cluster.
This is an opportunity for prior knowledge about the mixture distribution (<ref>) to enter the clustering procedure:
the score reflects the user's belief about how the mixture is structured.
SCAMP's default preference function prizes candidate clusters with low-variance, symmetric coordinate distributions.
SCAMP uses a combination of the dip test and the sample L-moments of <cit.> to quantify this preference.
If a user is analyzing data where they expect the mixture to have a different composition, the preference function can be modified to better suit specific domains. 
For example, the preference function might be modified so that candidate clusters with skew distributions are preferred.

Once the candidate clusters are scored, selection is cast as an instance of the maximum-weight independent set problem.
SCAMP attempts to pick a collection of candidate clusters with maximal total preference score, subject to the constraint that no two selected clusters contain a common observation of the data matrix.
SCAMP settles for a greedy solution to this problem (preferring an approximate-but-quick solution to refinements like those discussed in <cit.> and <cit.>).
This results in determining a subset of selected clusters from the set of candidate clusters.

This selection procedure often results in a large collection of residual observations that are not elements of any selected cluster.
From the perspective of the underlying data model, the residual observations are thought to be samples from different components of the mixture (<ref>).
To cluster these residual observations, the SCAMP search is repeated on the corresponding residual sub-matrix. 
This causes a new set of candidate clusters to be found, scored, and selected. 
These iterations continue until each row of the initial data matrix is an element of a unique selected cluster.

The final stage of the SCAMP procedure is annotation.
Since the default preference function causes SCAMP to prefer clusters with large dip-test p-values,
SCAMP assumes every selected α-m-cluster is unimodal along each of its coordinates.
SCAMP annotates a selected cluster by comparing the median value along each of its coordinates with a collection annotation boundaries for the same coordinate (see section <ref> for further details). 

When the scale of measurement of a coordinate is meaningful, the annotation procedure can create clusters with interpretable labels.
Even when the coordinates are not directly interpretable (as is the case when SCAMP is applied to principal components), this annotation stage has the effect of merging clusters.
This merging can be justified in view of the data model: merging two selected clusters indicates that they were two independent samples from the same component of the mixture (<ref>).
The clustering returned by a single SCAMP iteration is made up of a disjoint collection of labeled clusters.



 §.§ Outline of Paper


The remainder of this paper is organized as follows.
In section <ref>, we develop an extended version of the dip test of <cit.> that is used by SCAMP to search for candidate clusters containing a small number of observations.
Section <ref> provides the details of the SCAMP clustering algorithm.
The performance of SCAMP is examined in section <ref>, where we apply it to simulated and real data sets.
Proofs of several claims made in section <ref> are contained in the appendix.

As mentioned earlier, SCAMP often uses the taut string of <cit.> to induce modal groups.
We decided to use the taut string since, as noted by <cit.> page 1099, observation (iii), it has the minimum modality over a large collection of density estimators.
We found this modal sparsity appealing given SCAMP's search strategy.
However, SCAMP always uses the default value of κ=19 (discussed and set by the authors in <cit.>) to solve the κ-Kuiper problem.

For data matrices with fewer than 400 observations, 
the splitting is done with an extended version of the dip test in a sequential procedure designed to estimate the number of modes k present in the coordinate of the sub-matrix under consideration. 
In part, our extension of the dip test grew out of a desire to minimize SCAMP's dependence on bandwidth selection.
Once estimated, we use k̂ in conjunction with the one-dimensional dynamic programming version of k-medoids by <cit.>
to induce modal sub-collections (and so ensure the exactly k̂ modal sub-collections are produced). 

To our knowledge, this extended version of the dip test is new.
Other testing procedures, such as the excess mass test, could have been used in its place:
<cit.> discuss additional testing procedures.
However, as noted by <cit.> as well as <cit.>, the excess mass test concludes the presence of two modes in situations like those depicted by figure <ref>.
SCAMP is designed with the goal of detecting such small groups and selecting them if they have high preference scores: using our extension makes this outcome possible.
We also considered using the multiscale procedure of 
<cit.>
and the jaws-based method of 
<cit.> to test for multimodality in small samples.
We decided against using these tests since they both required developing automated procedures to make decisions on the basis of interval collections:
in the former, the collections 𝒟^+(α) and 𝒟^-(α);
in the latter, W- an M-components contained within larger candidate antimodal sections of the empirical distribution.
Thinking about how to develop such automated procedures led us to consider our extension to the dip, which we discuss in the next section.





 §.§ Summary of Contributions


In this paper, we have developed a new clustering algorithm called selective clustering annotation using modes of projections (SCAMP). SCAMP is designed to be:


  
  * Interpretable. Clusters produced by SCAMP are automatically assigned labels in terms of the columns of the underlying dataset.
      When the measurements in the underlying dataset are themselves meaningful, this produces clusters that human users can understand simply by reading their label.
  
  * Robust. SCAMP defines clusters in terms of the empirical distributions of their coordinate projections: they must be unimodal.
      SCAMP is consequently able to produce clusters with long tails that accommodate outliers.
  
  * Heterogeneous. SCAMP is able to produce clusterings of data matrix consisting of clusters of different sizes.
      
      
      
  
  * Simple to use. SCAMP determines the number of clusters in a dataset  as a consequence of three tuning parameters: a level α, 
      against which p-values from dip test are compared; a lower-bound m,
      which defines the smallest number of observations that can constitute a cluster; and a 
      univariate Gaussian variance parameter γ.
      The parameter α describes a univariate quantity, no matter the dimension 
      of the data matrix.
      Given these parameters, SCAMP will cluster the dataset:  
      estimating the number of clusters in a data matrix is an implicit part of its clustering 
      strategy.
      This is in contrast to procedures such as k-means and k-medoids, 
      which either require the user to set the number of clusters k themselves, or to estimate 
      this parameter prior to clustering.
  
  * Composable. The most computationally demanding task in the SCAMP procedure is the 
  search for candidate clusters.
      Each tree grown in the search requires computing the dip test and taut string density 
      estimator numerous times. 
      Fortunately, O(n) algorithms exist for both computations.
      In addition, the search can be parallelized: a parallel   implementation of this 
      procedure with R interface is available.
      When combined with a random-sampling scheme, it is possible to use SCAMP to cluster datasets 
      with a large number of observations
      and moderate number of features. We demonstrate this possibility on both real and simulated 
      data in section <ref>.
      In modern datasets, where the number of features is much larger than the number of 
      observations, 
      SCAMP can be used productively 
      in conjunction with dimensionality reduction methods such as 
      PCA.
      We demonstrate this possibility on bulk RNA seq data in section <ref>.


 



§ DOUBLE DIPPING


To introduce our extension of the dip test, we start with an informal discussion of the dip itself.
The basic observation underpinning the dip test of <cit.> is that the distribution function (df)
of a continuous unimodal distribution (suitably defined) has an s-shape.
Suppose a unimodal df with support on ℝ has a unique mode m.
On the interval (-∞,m]  such a df is convex; on [m,∞), it is concave. Taken together, they form the s-shape. 

Independent samples from a unimodal df will produce empirical dfs that resemble the underlying s-shape with increasing fidelity as the sample size grows larger.
The dip takes advantage of the fact that the greatest convex minorant (GCM) and least concave majorant (LCM) of an empirical distribution function:

  
  * Are inexpensive to compute.
  
  * Do not require the user to select any tuning parameters in order to compute them.
  
  * Can be combined with a linear component in order to estimate a sample's departure from unimodality under the null hypothesis that they originate from a unimodal (s-shaped) df.


Now, suppose we are examining a bimodal df F with support on ℝ that has two unique modes m_1, m_2, and antimode c_1, such that m_1 < c_1 < m_2.
The idea behind our extension is the observation that such a distribution function consists of two unimodal sub-dfs stitched together:
on the interval (-∞,m_1] it is convex; on [m_1,c_1], concave; on [c_1,m_2] it is again convex; and concave again on [m_2,∞).
Restricting to the interval (-∞,c_1] yields one s-shaped sub-df; restricting to [c_1,∞) another.
See figure <ref> for an example.



Continuing with this bimodal F, we define

    H_1 ≡  H_(-∞,c_1]= F/F(c_1)· 1_(-∞,c_1)+1_[c_1,∞)

and

    H_2 ≡ H_[c_1,∞)=F-F(c_1)/1-F(c_1)· 1_[c_1,∞) .

Both H_1 and H_2 are unimodal dfs in their own right.
From this it follows that the 1-dip will tend to 0 if applied to data from either H_(-∞,c_1] or H_(c_1,∞).
This suggests the 1-dip algorithm of <cit.> might be used to compute the 2-dip as well: we expect 
the maximum estimated dip of samples from the restrictions H_(-∞,c_1] and H_(c_1,∞) to be small.
Notice that this observation extends to an arbitrary (fixed) number of modes: any multimodal distribution (suitably defined) can be split at its antimodes into unimodal component dfs (when scaled).
Moreover, the dip can be applied to each of these restrictions in turn.
This suggests the following extension of the dip test.
Under the hypothesis of data being generated from a distribution with k-modes, split a sorted sample into k sub-collections at its k-1 antimodes and apply the 1-dip to each sub-collection.
If the k-mode hypothesis is true, we expect each of these restricted dips to be small. 

The remainder of this section extends the theorems of 
<cit.> to justify the preceding informal development. We extend the theorems under the assumption the location of antimodes are known.
Under such an assumption, the theorems of <cit.> extend naturally: we essentially inherit the original arguments.
In fact, most of the work is notational, in order to accommodate a distribution with multiple modes.
We develop necessary notation for the remainder of this section (taking the notation directly from <cit.>),
and conclude with a heuristic development of the algorithm we use to estimate the N-dip in practice.

For bounded functions F, G on ℝ define ρ(F,G) ≡sup_x |F(x)-G(x)|. For a collection of bounded functions 𝒜, define

    ρ(F,𝒜) ≡inf_A ∈𝒜ρ(F,A) .

Two classes have special importance.


  For N ∈ℕ, let 𝒰_N be the set of all distribution functions F on ℝ
  where there exist real numbers
  
    m_1 ≤ c_1 ≤ m_2 ≤…≤ m_N-1≤ c_N-1≤ m_N

  such that F is convex on (-∞, m_1], concave on [m_N, ∞), convex on [c_j-1,m_j] with
  2 ≤ j ≤ N, and concave on [m_j,c_j] with 1 ≤ j ≤ N-1.



  For N ∈ℕ, let 𝒱_N be the set of all functions f on ℝ
  where there exist real numbers
  
    0≡ c_0 ≤  m_1 ≤ c_1 ≤ m_2 ≤…≤ m_N-1≤ c_N-1≤ m_N ≤ c_N≡ 1

  such that f is a finite constant on (-∞, 0], (a possibly different) constant on [1, ∞),
  non-decreasing on (0,1),
  convex on [c_j-1,m_j] with 2 ≤ j ≤ N and concave on [m_j,c_j] with 1 ≤ j ≤ N-1.


We will extend the dip relative to the first class 𝒰_N.


  For fixed N ∈ℕ, define the N-dip of a distribution function F to be
  
    D_N(F) ≡ρ(F,𝒰_N) .



Note that when N=1, the preceding definitions reduce to those in the paper of <cit.>.
Also note the classes are nested: for all N ∈ℕ, 
𝒰_N ⊂𝒰_N+1
and 𝒱_N ⊂𝒱_N+1.
For the remainder of this section, suppose N ≥ 2, N∈ℕ, and that N is fixed.
The observations made in <cit.> persist for the N-dip: for distributions F_1,F_2 we have

    D_N(F_1)
      
      = inf_G ∈ 𝒰_𝒩sup_x |F_1-F_2+F_2-G| 
      ≤ D_N(F_2) + ρ(F_1,F_2) .

Moreover, for F ∈𝒰_N we have D_N(F) = 0, so D_N(F) measures the departure of F from distributions with N modes.
In what follows, let 𝒞^+_[a,b] denote the collection of all convex functions on the interval [a,b], and
𝒞^-_[a,b] the concave functions on the interval [a,b], with a,b ∈ℝ, a<b.
If either a= -∞ or b=∞, open the interval to extend the notation. We will use the GCM and LCM to measure departures from a modal hypothesis.


  The greatest convex minorant (GCM) of a function f on an interval [a,b] is defined
  
    GCM(F,[a,b]) ≡sup{G(x) | G(x) ≤ F(x)  on  [a,b],  and  G(x) ∈𝒞^+_[a,b]} .




  The least concave majorant (LCM) of a function f on an interval [a,b] is defined
  
    LCM(F,[a,b]) ≡inf{G(x) | G(x) ≥ F(x)  on  [a,b],  and  G(x) ∈𝒞^-_[a,b]} .



With these definitions in place, we first state a bound on the N-dip.


  For N ≥ 2, N ∈ℕ we have for any distribution function F on ℝ, we have
  
    D_N(F) ≤1/2N .



Proof of Lemma <ref>, along with the remaining claims in this section, are given in Appendix <ref>.
The proofs themselves are natural extensions of those given by <cit.>.

Our development and presentation of these extensions were also shaped by unpublished course notes of <cit.>,
available on his website.


We next observe for distributions concentrated on [0,1], departures from
the classes 𝒰_N and 𝒱_N are equivalent.


  Let N ∈ℕ be fixed. Let F be a distribution function with F(0)=0 and F(1)=1. Then
  
    ρ(F,𝒰_N) = ρ(F,𝒱_N) .



As with the dip, we next see that scaling and mixing the Uniform distribution with suitably bounded function scales the N-dip.


  Let F be a bounded function constant on [-∞,0] and on [1,∞]. Let I be the distribution
  function of the Uniform on (0,1). Let N ∈ℕ. Then
  
    D_N(α F + β I) = α D_N(F)

  for α, β≥ 0. 


The first theorem states that, as with the dip, the Uniform distribution is also extreme for the N-dip.


  Let N ∈ℕ.
  Let F_n be the empirical distribution function for a sample of size n from the Uniform on (0,1), and let B be
  the Brownian Bridge process with cov[B(s),B(t)] = s(1-t) for 0 ≤ s ≤ t ≤ 1. Assume that B is zero
  outside (0,1).  Then
  
    √(n)D_N(F_n) →_d D_N(B)  as  n ↗∞ .



The next theorem states that the N-dip, when applied to samples from multimodal distributions with exponential changes between modes and exponential decrease away from extreme modes, will tend to zero in probability.



  Let N > 2, N ∈ℕ be fixed. Let F ∈𝒰_N, and have non-zero
  k^th derivative at the distinct antimodes c_1,…,c_N-1 and the
  distinct modes m_1,…,m_N, for some k ≥ 2. For each ϵ, ϵ_1,
  ϵ_2 > 0, and for j ∈{1,2,…,N} that matches 
  the number of modes and antimodes, suppose
  
    inf_x < m_1-ϵ    d/dxlog F'(x)  > 0  ,  
    inf_x > m_N+ϵ    d/dxlog F'(x)  > 0  ,  
    inf_c_j-1+ϵ_1 < x < m_j-ϵ_2    d/dxlog F'(x)  > 0  ,  and 
    inf_m_j+ϵ_1 < x < c_j+1-ϵ_2    d/dxlog F'(x)  > 0  .

  Then √(n)D_N(F_n) →_p 0.


To summarize our findings, theorem <ref> suggests the N-dip will be 
able to detect a large number of distributions asymptotically.
Theorem  <ref> suggests the Uniform distribution may serve well as the null
distribution in finite samples. 
How do we compute the N-dip in finite samples? Theorem 6 of <cit.> 
and the subsequent discussion provide an algorithm for computing the 1-dip
[<cit.>].
Here we suggest an extension to estimate the N-dip. For notational clarity, we focus
our discussion on the 2-dip.

Returning to the bimodal example that began the section, we again consider the
restrictions (<ref>) and (<ref>).
Observe that

    F = H_1· F(c_1) +H_2· (1-F(c_1))  .

If we knew the location of the antimode c_1, the empirical distribution could be similarly partitioned:

    F_n = H_n,1· F(c_1) +H_n,2· (1-F(c_1))  ,

with H_n,1 and H_n,2 formed by F_n replacing F in the definitions of H_1 and H_2.
But, subject to the conditions of  theorem <ref>, √(n)D_1(H_n,1) →_p 0 and √(n)D_1(H_n,2) →_p 0.
Combined with Lemma <ref>, we infer

    √(n)D_2(F_n) 
       = √(n)· F(c_1)· D_2(H_n,1)  ∨√(n)· (1-F(c_1))· D_2(H_n,2)  
       = √(n)· F(c_1)· D_1(H_n,1)  ∨√(n)· (1-F(c_1))· D_1(H_n,2)  →_p 0 ,

where the second equality is suggested by the fact that
𝒱_1 ⊂𝒱_2.
Note the  order statistics X_(1),X_(2),…,X_(n) are consistent for their population counterparts: as the sample size increases a sample point will grow arbitrarily close to c_1.
In finite samples, we lack knowledge of F and c_1. Under the hypothesis of bimodality, plugging in empirical estimates leads naturally to 

    √(n)D̂_2(F_n) = √(n)· F_n(ĉ_1)· D_1(H_n,1)  ∨√(n)· (1-F_n(ĉ_1))· D_1(H_n,2) .

Since we do not make assumptions about the relative magnitude or location of the modes, the estimate ĉ_1 can be obtained by exhaustive data splitting.
Hence, for j ∈{4,5,…,n-3},  let ĉ_1,j = X_(j). 
The hypothesis of bimodality indicates that splitting the sample at the order statistic closest to the antimode will lead to scaled unimodal distributions that each minimize the maximum 1-dip of H_n,1 and H_n,2.
Our proposed estimator of the 2-dip is consequently

    √(n)D̂_2(F_n) = min_j ∈{4,5,…,n-3}[√(n)· F_n(ĉ_1,j)· D_1(H_n,1)  ∨√(n)· (1-F_n(ĉ_1,j))· D_1(H_n,2)] .

The restriction on j follows because the 1-dip achieves its lower bound with high probability on small sample sizes 4 (see Proposition 1 of <cit.>).
Conveniently, this formulation allows us to use the 1-dip algorithm of Hartigan and Hartigan to estimate the 2-dip. 
The formulation also extends naturally to the estimation of the N-dip, by partitioning order statistics of an observed sample into N-components (with each component containing at least 4 observations),
weighting them according to their empirical mass, and then searching for the sample index that minimizes the maximum over the collection.

For fixed sample size, theorem <ref> suggests the Uniform can be used to obtain p-value.
We use the Uniform distribution in practice: we have computed the distribution of
values obtained by the estimator on random samples from the Uniform distribution for
sample sizes up to N = 400 for the 2-dip and 3-dip. 
These values, along with the implementation of the 2-dip and 3-dip estimators, are
available in the package scamp.
Pseudo-code describing the computation of the 2-dip-estimate is given in figure (<ref>). 





§ SELECTIVE CLUSTERING ANNOTATED USING MODES OF PROJECTIONS
 



SCAMP is an iterative procedure: each iteration determines a labeled partitional clustering of a data matrix X.
To understand its iterative nature, we must first describe details of a single SCAMP iteration that were omitted in section <ref>.
We do so now.



 §.§ Single Iteration


Before SCAMP is applied to a data matrix X, there is a pre-processing step: the data matrix is normalized so each column vector has mean zero and unit variance.
After normalization, the data matrix is passed to SCAMP.

There are four phases in a single SCAMP iteration. 
In the first phase, two forms of structured noise are added to the data matrix.
This is done so that ties within column vectors of a data matrix are broken and the relative order of observations within a column are perturbed.
In the second phase, the data matrix is recursively searched for α-m-clusters.
In the third phase, a partial clustering of the data matrix is determined by selecting α-m-clusters with high preference scores.
If any observations are not assigned to a cluster after phase three,
the second and third phase repeat on residual sub-matrices until a complete partitional clustering of the data matrix has been found.
In the fourth and final phase, the selected clusters are assigned descriptive labels.
This labeling provides a name for every row of the data matrix.
Figure <ref> visualizes these four phases, when SCAMP is applied to the famous iris data set of <cit.>. 



We now look at the details of each phase.



  §.§.§ Add noise to the data matrix


For 1 ≤ i ≤ k, denote the rows of a data matrix X_k × p by 

    r_i ≡ (x_(i,1),…,x_(i,k)) ,

for 1 ≤ j ≤ p, the columns by

    c_j ≡ (x_(1,j),…,x_(k,j))'  ,

and specific entries as x_(i,j). 
A single scamp iteration begins by adding noise to each column c_j of the data matrix.
The noise procedure is the same for each column c_j. 
First, the order statistics of c_j 

    O_(1:d_1), O_(2:d_2),…,O_(q:d_q)

are determined. 
Here O_(1:d_1) indicates the minimum value observed in column c_j with d_1 repetitions, O_(q:d_q) the maximum with d_q repetitions, and O_(n:d_n) 
for 1 < n < q  the values between. Once determined, duplicate entries of c_j are replaced by samples from a uniform distribution.
Suppose an entry x_(i,j) corresponds to the order statistic O_(n:d_n). If d_n=1, x_(i,j) is unmodified. Otherwise x_(i,j) is replaced by

    u_(i,j)∼Uniform(O_(n:d_n)+O_(n-1:d_n-1)/2,O_(n:d_n)+O_(n+1:d_n+1)/2) .

Notice that d_n entries of c_j are replaced by samples from this specific uniform distribution.
In the case that either the minimum value of c_j is not unique, the lower bound of (<ref>) is replaced by that minimum value.
The same hold if the maximum value of c_j has copies, though with the upper bound of (<ref>) replaced by the maximum value instead.

The choice of adding uniform noise is motivated by SCAMP's reliance on the dip test. 
As discussed in section <ref>, the dip test is calibrated against the uniform distribution.
By replacing tied observations with draws from a uniform distribution that is supported on an interval bounded by neighboring order statistics, 
we have replaced a jump of d_n/k in c_j's empirical distribution with d_n jumps of 1/k. 
From the point of view of the dip test, the empirical distribution of c_j now behaves like a unimodal distribution on the interval 
([O_(n:d_n)+O_(n-1:d_n-1)]/2,[O_(n:d_n)+O_(n+1:d_n+1)]/2). 

Once uniform noise has been added to c_j, the order statistics 

    N_(1:1), N_(2:1),…,N_(k:1)

are again determined. Notice each order statistic is now unique. 
The noise-step concludes by associating each entry 
u_(i,j) with its new order statistic N_(n:1), and replacing it by

    n_(i,j)∼Normal(μ = u_(m,j) , σ = N_n+1:1+N_n-1:1/2·γ) ,

with γ a tuning parameter (by default we set γ≡ 4).
As almost all the mass is concentrated within three standard deviations of the mean, the addition of Gaussian noise (<ref>) provides
observations with the chance to swap relative position with some of their neighbors. 
For observations with values near modes of c_j, this step is largely inconsequential in their ultimate cluster assignment. 
But for observations near antimodes of c_j, this step increases the chance they change cluster assignment on successive SCAMP iterations.
Over multiple SCAMP iterations, changes in cluster assignments provide an uncertainty measure of each observation's SCAMP label.



  §.§.§ Search for candidate clusters


After structured noise has been added to the data matrix, SCAMP searches the data matrix for subsets of rows that follow the definition of an α-m-cluster.
To conduct this search, two parameters must be set: the significance level for the dip test of  <cit.>, α, and the lower bound on cluster size, m. 

Once set, SCAMP begins by testing each coordinate projection of the dataset for multimodality using the dip test.
For each column in the data matrix, when multimodality is detected by the dip test
modal groups are induced by estimating the density of the column vector and then separating observations into groups relative to the antimodes of the density.
This is done using the taut string density estimator of <cit.>.
The cut-points that separate the modal groups are determined by calculating the average coordinate value along each antimodal component of the taut string.
Sub-matrices are then created for all observations that fall between each neighboring pair of cut-points (the minimum and maximum value along the coordinate sub-collection are treated as cut-points).
Within each sub-matrix, projections of the remaining p-1 coordinates are recursively tested for multimodality (again using the dip test at level α).

If multimodality is detected in a coordinate of a sub-matrix containing fewer than 400 observations, the taut-string is not used to induce modal groups.
Instead, bimodality is tested using the 2-dip of section <ref> at level-α/2.
Failure to reject the hypothesis of bimodality causes SCAMP to conclude k̂=2 modal sub-collections are present in the coordinate.
If bimodality is rejected, trimodality is tested using the 3-dip at level-α/3.
Now, a failure to reject leads SCAMP to conclude k̂=3 modal sub-collections are present.
In principle, this sequential testing procedure can be carried out until it arrives at an estimate of n-modality.
However if tri-modality is rejected, SCAMP simply concludes k̂=4 modal sub-matrices are present due to computational limitations.
SCAMP then separates the projection into modal sub-collections using
the one-dimensional version of k-medoids by <cit.> with k̂ groups.

We introduce this secondary sequential testing procedure in sub-matrices with fewer than 400 observations for several reasons.  The first is computational. As the number of observations increases, the cost of exhaustive splitting grows with the falling-factorial of the number of modes of the N-dip. Because of this computational cost, we have not yet been able to compute tables of critical values for the N-dip when N > 3. 

Ideally, we would like to apply this sequential procedure until the N-dip fails to reject its null hypothesis for some N, in place of the current combination of the dip test and subsequent density estimate. This is because the sequential procedure relies only on the choice of α. It does not introduce dependence on bandwidth into the SCAMP search procedure, implicitly or explicitly. The taut string, for example, requires setting the κ parameter, or accepting its default value of κ = 19. 

In our experience, when applying SCAMP to real datasets the default κ-value of 19 causes the taut string to produce very useful density estimates the vast majority of the time. However, in rare data sets we have found that the taut string can find a very large number of antimodes when applied to certain features. In these cases, exhaustive search of the annotation forest becomes impossible unless the feature is modified or removed from the analysis. 

To account for this rare circumstance, our implementation in the scamp package has a parameter to
limit the number of antimodes SCAMP will search along any coordinate in a given partition tree. In the 
event the taut string estimate exceeds the parameter value, the search terminates along that branch. In
our typical use cases, we set this parameter to 100, though we expect this parameter should be
modified depending on the problem domain.

If we view an estimated density having a very large number of antimodes when in truth there are very few antimodes in the underlying distribution as an instance of type I error, we can think of switching to the sequential testing procedure as explicitly imposing a preference for type II error into the SCAMP search procedure in small sub-collections. The current implementation guarantees that SCAMP will never produce greater than 4 modal sub-collections when examining small subsets of the data matrix, even in cases where there truly are more than 4. Should more efficient implementations of the N-dip become available, this restriction can be reduced for larger number of modes in the N-dip and more than 400 observations in a sub-matrix.

The recursive search terminates along a branch when it encounters a sub-matrix whose coordinates all have dip test p-values exceeding α,
or the sub-matrix under consideration contains fewer than m observations. 
In the former case, SCAMP has found a candidate cluster and the indices of the rows in the sub-matrix are recorded.
In the latter case, stopping is viewed as uninformative and the indices are deleted.

Ideally this search step is exhaustive. 
In practice, exhaustive search is impractical for data matrices with even a moderate number of
features.
Sampling is used to address this limitation: our implementation in the package scamp is able
to search uniformly at random among the space of partitional trees in the annotation forest until
finding a specified number of candidate clusters. We empirically justify this sampling approach with
our simulated data experiment in section <ref>.



  §.§.§ Select Candidate Clusters


Clustering of the dataset begins with SCAMP imposing a preference values on the set of 
candidate clusters collected on phase 2.
As discussed in section <ref>, the preference function reflects the users beliefs 
about components in the mixture (<ref>).
For its default, SCAMP prefers candidate clusters that are unimodal, symmetric, and have low variance 
along each coordinate projection.
If SCAMP is applied to a data matrix produced in a domain where other cluster shapes are desirable, the
preference function can be modified to select accordingly.

SCAMP uses 
the sample L-moments of <cit.> to assign preference
scores to candidate clusters.
Recalling definition <ref> from section <ref>, 
denote the set of s candidate clusters found in SCAMP phase 2 by

    𝒞≡{X_ℐ× p^1,…,X_ℐ× p^r,…,X_ℐ× p^s } ,

with 1 ≤ r ≤ s and with the indices into the original data matrix different for each element of 𝒞.
Denote Hosking's trimmed second, third, and fourth order
for the i^th  coordinate of X_ℐ× p^r by σ_i^r, τ_i^r, and ϕ_i^r  respectively, with 1 ≤ i ≤ p. 
Additionally, denote the p-value of the dip statistic D̂_1  for each of the p-coordinates of of X_ℐ× p^r by δ_i^r.

For each fixed candidate cluster X_ℐ× p^r in 𝒞, SCAMP assigns it a preference value by first computing

    δ^r ≡min_1 ≤ i ≤ pδ_i^r  , σ_max^r ≡max_1 ≤ i ≤ pσ_i^r , σ_sum^r ≡∑_1 ≤ i ≤ pσ_i^r ,
       τ_max^r ≡max_1 ≤ i ≤ pτ_i^r ,  τ_sum^r ≡∑_1 ≤ i ≤ pτ_i^r , ϕ_max^r ≡max_1 ≤ i ≤ pϕ_i^r ,  and ϕ_sum^r ≡∑_1 ≤ i ≤ pϕ_i^r .

With these quantities computed for each element of 𝒞, SCAMP next computes

    σ_max≡max_1 ≤ r ≤ sσ_max^r   , σ_sum≡max_1 ≤ r ≤ sσ_sum^r  ,
    τ_max≡max_1 ≤ r ≤ sτ_max^r   , 
        τ_sum≡max_1 ≤ r ≤ sτ_sum^r  ,
    ϕ_max≡max_1 ≤ r ≤ sϕ_max^r   ,  and  ϕ_sum≡max_1 ≤ r ≤ sϕ_sum^r  .

Using these maximums across 𝒞, normalized values of the sample L-moments are then computed for each X_ℐ× p^r: 

    σ_nm^r ≡ 1 - σ_max^r/σ_max , σ_ns^r ≡ 1 - σ_sum^r/σ_sum ,

along with the analogous quantities τ_nm^t, τ_ns^r, ϕ_nm^r, and  ϕ_ns^r.
The preference score of cluster X_ℐ× p^r is then defined to be

    𝒫(X_ℐ× p^r) ≡δ_r + σ_nm^r+σ_ns^r/2+ τ_nm^r+τ_ns^r/2+ ϕ_nm^r+ϕ_ns^r/2 .

The normalization step (<ref>) allows the following interpretation of the preference score (<ref>): the higher the better. 
Notice that X_ℐ× p^r enters 𝒞  because δ_r > α. 
By using δ_r as a signal of unimodality in (<ref>), we  SCAMP is encouraged to discriminate the degree of unimodality.
Suppose SCAMP is applied to a data matrix with α = 0.05. 
Two hypothetical candidate clusters X_ℐ× p^r_1 and X_ℐ× p^r_2 enter 𝒞 with d_r_1 = 0.06 and d_r_2 = 0.94. 𝒫 
on the basis of the dip alone, 𝒫 is defined to prefer X_ℐ× p^r_2.
The preference score (<ref>)  averages the trimmed sample L-moments to balance between the worst shape and average shape of the coordinates of a candidate cluster X_ℐ× p^r. 
For example, (σ_m^r+ σ_s^r)/2 is meant to score a candidate cluster so that if either one single coordinate has extremely high spread or if many coordinates have moderate spread,
then SCAMP will not find the candidate cluster particularly appealing. Similar reasoning motivated (τ_m^r+τ_s^r)/2 and (ϕ_m^r+ϕ_s^r)/2.
The sample L-moments are trimmed in order to moderate the preference score penalty incurred by a candidate cluster when it contains an outlying observation along any of its coordinates.

By scoring each candidate cluster, SCAMP has created a set of weights

    𝒫≡{𝒫(X_ℐ× p^1),
    …,
    𝒫(X_ℐ× p^r),
    𝒫(…,X_ℐ× p^s) }

associated with the set of candidate clusters 𝒞. A graph can now be induced by viewing each
element of 𝒞 as a node, and inducing the adjacency matrix with edge set

    ℰ≡[ 
    𝔼(X_ℐ× p^r_1,X_ℐ× p^r_2)
      𝔼(X_ℐ× p^r_1,X_ℐ× p^r_2) = 
    
    0  if   X_ℐ× p^r_1∩ X_ℐ× p^r_2 = ∅
    
    1  otherwise
    ] .

That is, two candidate clusters are connected by an edge if they share a common observation from the data matrix. 

Clustering of the data matrix proceeds by selecting a subset of 𝒞 that maximizes the sum of the associated preferences 𝒫  and  no two selected clusters are connected in ℰ.
This is an instance of maximum weight independent set (MWIS) integer program.
The MWIS problem is NP-hard; guarantees on the performance of approximate solutions are not known expect in special cases [<cit.>].
SCAMP performs a greedy search: it sequentially picks the available candidate cluster with the highest preference score, 
and eliminates all candidate clusters connected in ℰ from contention in subsequent selection steps. 



  §.§.§ Cluster residual observations







The initial set of selected clusters, those candidate clusters picked by the greedy selection procedure of phase 3, often determine a partial clustering of the data matrix.
This is because the set of selected clusters consists of leaves of different partition trees in the annotation forest.

The residual observations that are not elements of any selected clusters do not necessarily obey the definition of an α-m-cluster.
When these residual observations exhibit multimodality at level α along some of their coordinate projections, they are thought to represent samples from components of the mixture (<ref>) that were not found in the initial search.
Therefore, a SCAMP iteration does not terminate after a single search-and-selection pass.

To complete the clustering of the data matrix, SCAMP conducts the search phase 2 and the selection phase 3 on the subset of rows of the data matrix corresponding to the residual observations.
This creates a new set of candidate clusters 𝒞. The new set is scored according to the preference function 𝒫. 
New clusters are then selected and appended to the clusters already selected in the earlier selection round.
Once again, the secondary selection procedure may produce residual observation.
So, the procedure recurses, applying the search and selection phases to smaller and smaller subsets of the original data matrix until either no residual observations appear, 
the residual observations comprise an α-m-cluster, 
or there are fewer than m residual observations.

At the completion of this process, a set of selected clusters 𝒮 has been created. 
By its construction, 𝒮 partitions the data matrix into α-m-clusters. 
This is why SCAMP does not require the user to select the number of clusters as a tuning parameter: the number of clusters follows instead from the choice of α and m. 

α may initially seem as difficult a parameter to set as the number of clusters k.
In practice, however, we have found it relatively simple to set since α describes
a one-dimensional characteristic of the data matrix.
Our standard procedure is the following.
Before we apply SCAMP to a new data matrix for a given application, we first
decide on an absolute upper bound on α. 
Intuitively, this corresponds to the minimum amount of evidence we require to split 
a coordinate into modal sub-collections.
Based on our experience of applying the dip test to a wide-range of datasets, we set this 
upper bound to a default value of α = 0.25.

Next, we refine this upper bound on α by  
adding the structured noise of phase 1 to the columns of the data matrix, then computing the 
dip test p-values for each column vector, and finally plotting their distribution.
If there is a clear point of separation in the distribution of p-values that
is slightly larger than our default choice of α=0.25, 
we can modify our choice upwards to respect the observed separation.
Our case studies in section <ref> provide several examples 
of these plots and how we use them.

When there is no clear point of separation in the p-values but a subset of the p-values fall below
our default α value, we normally use the default.
However, if we want a sparser clustering of the data matrix, we 
next look at histograms of the features below our bound,
and pick an α < 0.25 value large enough to ensure those features which appear 
multimodal from inspection are entered into the search phase.

In our experience, classical default values, such as α = 0.05, can also work well for
high-dimensional datasets with a large number of observations, 
so long as the signal of interest is well-separated:
the separated mixture simulation of section <ref> is one example of 
such a scenario. 
To summarize: while the choice of α can vary from problem to problem, 
inspecting the distribution of dip-test p-values for the available features provides a simple
data-driven approach for modifying this parameter from a default setting of α = 0.25.

As for the parameter m, if prior knowledge informs the expected cluster size, we can use this
knowledge to set a value. If no prior knowledge informs the analysis, we simply pick as small a
value as computational resources allow: our default value is m=25.
Beyond the computational limit, there is an absolute lower bound of m = 4 due to SCAMP's
reliance on the dip test.



  §.§.§ Annotate selected clusters


A SCAMP iteration concludes by assigning labels to the selected clusters 𝒮. 
Since each selected cluster is an α-m-cluster, we expect the distribution of observations within a cluster to be somewhat unimodal along any given coordinate.
Empirically we have observed that while the hypothesis of unimodality generally holds for clusters selected early in the selection process,
clusters selected near the end of the selection process can be asymmetric with shoulders.
Nevertheless, our annotation procedure behaves as if the assumption of unimodality holds equally well for all selected clusters.

To begin, annotation boundaries are determined for each of the coordinates of the data matrix X_k × p that are multimodal at level α.
The annotation boundaries are computed using data collected in the annotation forest search of phase 2.  
For each of the j coordinates of the data matrix, 1 ≤ j ≤ p, SCAMP records the location of all cut-points found in all trees of the annotation forest grown during the initial
candidate cluster search. It records each collection of cut-points in a separate set. Define

    K_i^j ≡{(κ_1,…,κ_i)   SCAMP induces  i+1 modal groups in coordinate  j
      with i cut points κ_1 < … < κ_i
    } ,

and the cardinality of each set by 𝒦_i^j = K_i^j. Define m^j  to be the largest index  i  such that 𝒦_i^j > 0.  
The annotation set for the j^th coordinate is K_a^j, with a defined to be the smallest index 1 ≤ a ≤ m^j such that 𝒦_a^j  ≥𝒦_k^j for all 1 ≤ k ≤ m^j.

Supposing the annotation set K_a^j has cardinality 𝒦_a^j = n,
the annotation boundaries K^j are derived for the j^th coordinate by computing the median cut-point 
for each of the a annotation coordinates:

    K^j ≡{(k_1,…,k_a)    k_b ≡median_1 ≤ i ≤ n κ_(b,i) for  1 ≤ b ≤ a } .

To annotate a specific coordinate in a selected cluster, three sample percentiles of the selected cluster are computed. 
By default the q_l ≡ 50^th - ϵ, q_m ≡ 50^th, 
and q_u ≡ 50^th + ϵ percentile are used, 
leading us to annotate clusters relative to their median coordinate values.
These values can be adjusted by the user, depending how rigidly they want the final labels 
assigned to a cluster to respect the annotation boundaries along all labeled coordinates.
In our experience, this varies from problem to problem.
The label provided to each cluster is a proportion describing the relative position of the sample 
quantiles to the annotation boundaries (<ref>),
thereby associating the cluster with a mode of the coordinate projection.

In the special case a=1, the following occurs: if q_l ≥ k_1, 
the cluster is given the label 1/1; if q_u ≤ k_1, the cluster is given the label 0/1;
otherwise, the cluster is split, with points in the cluster below k_1 given the label 0/1, 
and those above k_1 the label 1/1.
If a > 1, then q_m is compared to each of the boundaries k_1,…,k_a.
If q_m < k_1, the cluster is labeled 0/a; otherwise the cluster is labeled with the index of 
the largest annotation boundary q_m exceeds.

To summarize: for each multimodal coordinate of the data matrix X, SCAMP clusters are assigned 
labels with fractions describing their relative position of the cluster column 
sample quantiles to annotation boundaries. 
In the scamp package, 
these fractions (up to 8) are then mapped to a dictionary to give each 
observation an interpretable label for the particular fraction. 
For example, if a=1, SCAMP reports 0/1 as “lowest”, and 1/1 as “highest”; if a=2, 
1/2 is “medium”; etc.  



 §.§ Multiple Iterations


As noted at the start of this section, a single SCAMP iteration determines a partitional clustering of a data matrix X_k × p.
The clustering produced by that iteration is always dependent on the realization of the random noise in section <ref>.
Further uncertainty is added if the search of section <ref> finds only a random subset of the annotation forest.

SCAMP is able to quantify the effect of these stochastic sections of the procedure by being run
multiple times on the same data matrix X_k × p.
By providing a name for each selected cluster in a given iteration the annotation step, SCAMP also
provides names for each individual observation in the data matrix.
Because of this, we can track the naming history of each observation in the data matrix across
multiple iterations.
We emphasize the computation efficiency: 
keeping this record only requires updating a map for each observation across the iterations, 
rather than recording which observations are co-clustered as SCAMP iterates.
Labels assigned to a given observation define the keys for each map; the values are integer counts 
of the number of times the observations is labeled a particular way.

We view the name determined by a single SCAMP iteration  
as a vote for how a particular row of a data matrix should be described.
For some rows, the vote is consistent: every iteration SCAMP gives a specific
observation the same name. 
In such a case, we have high confidence the observation is named appropriately.
For other rows, the vote is split, and two or three names appear with some regularity. 
In these cases, it seems likely the appropriate name might be some combination 
of those names that appear regularly.
And for still other rows, the vote is inconclusive, with no collection of names appearing to be 
more than fluctuations of the SCAMP procedure itself.

Because of this, the output of the scamp procedure in the package scamp is not one 
but two clusterings of a data matrix X_k × p.
The first is the clustering found by labeling each observation according to label it is given most
often across multiple observations.
The second is the clustering found by a run-off heuristic, which we now describe.

If, across multiple SCAMP iterations, an observation receives the same label in more than half of
them, the observation is named according to that label.
However, if an observation is given a label less than half but more than 30% of the time, and
the second most frequent label appears more than 20% of the time,
then the two labels are combined according to the following heuristic. If the most frequent label
describes a coordinate, but the second most frequent label does not, 
the final label is given the most frequent label's description. However, if the most frequent
label describes a coordinate by the ratio a/b  and the second most frequent label 
described it as c/d, that coordinate is then called (ad + bc)/(2bd). For example, if the most
frequent label is 0/1 (corresponding to the annotation “lowest”), and the second most frequent
label is 1/2 (corresponding to the annotation “medium”), the combined label would become
(0· 2 + 1· 1)/(2· 1 · 2) = 1/4 (corresponding to the annotation “medium-low”).

A similar combination heuristic is carried out in the case where the most frequent label for an observation occurs between 20% to 30% of the time, and the second and third 
most frequent more than 15%. 

Of course, other voting schemes could be considered in place of these heuristics, and they carry
interesting questions of their own: is any one voting scheme optimal, or does it depend on problem
domain? We note that instead of running SCAMP for a pre-specified number of iterations, SCAMP
could instead run until subsequent iterations produce 
no change (larger than ϵ) in the clustering determined by the selected voting scheme. 
In our own use of SCAMP, we most-often use the result of the maximum-vote heuristic after setting
the number of iterations to as large a value as computation time allows.  





§ SIMULATIONS AND EXAMPLES



In this section we apply SCAMP to a variety of data sources, 
both real and simulated.
We use the adjusted rand index [<cit.>] 
as one measure to quantify clustering performance, and refer to it by ARI.
We use the VI distance of <cit.> as another. 

Two independent clusterings of the same data matrix 
have an ARI value of 0 in expectation. Two
identical clusterings have an ARI value of 1.  
<cit.> notes that the
operative base-line of the ARI varies from 0.05 to 0.95 (pg. 876) and may be negative when
comparing certain clusterings (pg. 886). 
On the other hand, <cit.> proves that
the VI distance is a metric on the space of clusterings of a data matrix, ranges in value between
0 (when clusterings are identical) and log n (with n the number of observations), and can
be normalized to compare clusterings across data matrices with different number of observations. 
We report both values when analyzing real data matrices. In our simulation studies, we report only
the ARI to save space in graphical summaries.

Unless otherwise noted, SCAMP uses its default values α=0.25, m=25, and γ=4
in all case studies.




 §.§ Iris Data


We begin by showing how SCAMP clusters the iris data set of <cit.>. 
The data matrix contains measurements of sepal length, sepal width,
petal length, and petal width for 150 irises of three types: setosa, versicolor, and virginica.
There are multiple ties due to rounding: only 43 petal length measurements are unique, 
the maximum of the four variables. 
Because of the numerous ties SCAMP's noise phase has a substantial effect on the clustering.

Before running SCAMP, we first look at the distribution of dip test p-values for the data matrix.
Figure <ref> shows this distribution.
Based on this distribution, we see that choosing smaller values of α will not change
which columns are annotated: petal length and petal width have dip test p-values below 0.01
every iteration.
We note, however, that more conservative α values will lead
to a sparser annotation forest, since the same α value is used to split a coordinate across
all depths of a partition tree.
We proceed with our default value of α = 0.25. 



Figure <ref> summarizes the clusterings found across 500 SCAMP iterations with
with α=0.25, m=25, and the Gaussian noise parameter (<ref>) 
γ = 4. 
Each iteration, SCAMP randomly samples 200 candidate clusters from the annotation forest.
The cluster labels only describe an observations petal width and petal length, 
since sepal width and sepal length appear unimodal to SCAMP at level α=0.25. The labels 
determined by SCAMP appear useful in characterizing the differences between the three iris 
species.





 §.§ Olive Oil


We next apply SCAMP to the olive oil dataset of <cit.>.
This dataset has been previously analyzed by <cit.>, <cit.>, and
<cit.>.
We include this example to show how the descriptive labels produced by SCAMP can be useful for data
analysis, since the derived cluster labels explain differences between observations in different SCAMP clusters.

A description of these data is provided in the R package classifly [<cit.>],
from which we quote:

The olive oil data consists of the percentage composition of 8 fatty acids 
(palmitic, palmitoleic, stearic, oleic, linoleic, linolenic, arachidic, eicosenoic) 
found in the lipid fraction of 572 Italian olive oils. 
There are 9 collection areas, 4 from southern Italy (North and South Apulia, Calabria, Sicily), 
two from Sardinia (Inland and Coastal) and 3 from northern Italy (Umbria, East and West Liguria).


The 9 collection areas constitute 3 production regions in this dataset: 
southern Italy, Sardinia, and northern Italy.

Before running SCAMP, we again check our default choice of α against the
the data matrix: figure <ref> shows the distribution dip test 
p-values of the features in the data matrix.
We see the top 6 features always have a dip test p-value below our default α=0.25
aross 1000 noise iterations.
We observe, however, that the stearic fatty acid feature has a dip test p-value slightly above 
0.25 almost 35% of the noise iterations. 
This poses a problem for SCAMP's annotation phase, since it indicates the stearic fatty acid feature
will often drop out of the cluster annotation strings. 
This in turn will bias the maximum vote annotation heuristic.



Because of this, we decide to increase our default α: we 
run 5000 SCAMP iterations with α=0.30, m=25, and γ=4. 
Each iteration, SCAMP randomly samples 400 candidate clusters from the annotation forest.
The maximum-label clustering heuristic determines 11 clusters after 5000 iterations.
Of these 11 clusters, only 8 have more than 25 observations, our value for m.
We emphasize SCAMP can produce clusters with fewer than m elements since the final 
clustering is determined by the most frequently assigned
label for each observation across the 5000 SCAMP iterations.
The parameter m only restricts cluster size in 
the search for candidate clusters for a single iteration. 

Using the 9 Italian collection areas as the true cluster assignment,
the maximum-label
SCAMP clustering has an ARI of 0.734 and unadjusted VI distance of 1.377. 
Each cluster is labeled according to the relative quantities of the 8 
measured fatty acids.
For example, one cluster of 190 observations is labeled 
“palmitic highest, palmitoleic highest, stearic lowest, oleic lowest, linoleic highest, linolenic 
highest, arachidic highest, eicosenoic highest”. A second cluster of 57 observation is labeled
“palmitic highest, palmitoleic highest, stearic medium-high, oleic lowest, linoleic lowest, linolenic
highest, arachidic highest, eicosenoic highest”. A third cluster of 34 observations is
labeled “palmitic lowest, palmitoleic lowest, stearic medium-high, oleic highest, linoleic lowest, 
linolenic highest, arachidic highest, eicosenoic highest”.
These clusters broadly correspond to the southern Italy production region, which
contains the collection areas  South-Apulia (206 observations), North-Apulia (25 observations), Sicily (36 observations), and Calabria (56 observations).
We visualize these clusters in figure <ref> using t-SNE [<cit.>] to map the eight dimensions to two.

By coloring the t-SNE map according to the relative magnitude of each olive 
oil's fatty acid content, we see that the SCAMP labels reflect differences in the 
underlying fatty acid measurements. For example, olive oils in the cluster “palmitic 
highest, palmitoleic highest, stearic lowest, oleic lowest, linoleic highest, linolenic
highest, arachidic highest, eicosenoic highest” are mostly concentrated in a region of
the t-SNE map with high and low measured values of these fatty-acids, respectively (relative to the 
data set). A similar observation holds for the other clusters determined by SCAMP. 

From a data analysis standpoint, this allows an analyst to immediately understand why 
a given olive oil ends up in a particular SCAMP cluster. 
For example, within the South Italy production region, the three main 
SCAMP clusters all contain olive oils with
relatively high amounts of linolenic, arachidic, and eicosenoic fatty acids. 
The two SCAMP clusters that correspond (largely) to South-Apulia and Calabria also have
olive oils with 
high amounts of palmitic and palmitoleic acid, while the SCAMP cluster corresponding
(largely) to North-Apulia and Sicily has olive oils with
low amounts of these same fatty acids. On the other hand,
the olive oils in
the SCAMP clusters corresponding to
North-Apulia and Calabria have high amounts of stearic and 
linoleic acides, while the olive oils in the
South-Apulia/Sicily cluster has low amounts of these same fatty acids. 









 §.§ Separated mixtures


Here we conduct a simulation study to model the following setting: an experimenter has
collected n observations for p variables. The experimenter wishes to cluster the data 
in order to find sub-collections of rows that exhibit a common signal along some sub-collection 
of the p columns. 
To generate a matrix with such a property, we sample each sub-collection 
of rows from a multivariate distribution.
In the simplest case, the multivariate distribution is a Gaussian. 
Both the mean vector and covariance matrix differ from cluster to cluster.
In this simplest case, "no signal" corresponds to a mean vector entry 0 and an
activation "signal" corresponds to a mean vector entry larger than 0.

We examine three general scenarios.



  * Scenario 1: Moderate sample size, moderate number of clusters.

  * Scenario 2: Large sample size, moderate number of clusters.

  * Scenario 3: Large sample size, large number of clusters.


There are five binary parameters that can be used to modify a scenario.
Together they create 32 distinct simulation settings within a scenario.
For each setting, we run multiple iterations of the simulation.
We provide a graphical overview here to illustrate how the simulation parameters change the underlying 
data matrix in figure <ref>.
The simulation parameters, and the modifications they determine, are described in detail in appendix 
<ref>



In addition to the graphical depiction of the simulation parameters provided in figure <ref>,
we also provide t-SNE visualizations of the cluster structure produced in
data matrices under different parameter settings in figure
<ref>.
t-SNE is particularly well suited to visualizing these data since many of the clusters 
are sampled from multivariate Gaussian or multivariate t distributions.
The t-SNE maps show significant overlap when the observations in clusters
are taken through the maps described in the appendix <ref>.

In our experience, overlapping clusters are commonly observed in t-SNE maps produced 
from biological datasets.
Looking ahead, we see such overlap in our analysis of the GTEx data, visualized in figure <ref>.





  §.§.§ Scenario 1: Moderate sample size, moderate number of clusters


The graphical overview of figure <ref> is specialized to this scenario.
In it, the data matrix contains 3000 observations for all 32 parameter settings.
There are either 10 clusters in 20 dimensions, or 17 clusters in 40 dimensions.
Cluster sizes range between 125 and 1000 observations in the 10 cluster scenario.
The clustering ground truth is taken to be the label of the mixture component which produces an observation of the data matrix.

We sample data matrices 10 times under each simulation setting. 
Each iteration, we apply the following clustering methods. 
Notice that several methods require the user to set the number of clusters as a parameter.
In the current simulation scenario, we provide all methods except SCAMP the number of clusters.



  * (Oracle) Affinity Propagation (AP) [<cit.>]: number of clusters provided to the procedure “apclusterK” each iteration.
    Data matrix scaled to mean 0 and unit variance.
    We set the parameters “maxits=2000”,“convits=2000”, and “bimaxit=100”.

  * (Oracle) K-Means [<cit.>]: number of clusters set to the truth each iteration. Data matrix scaled to mean 0 and unit variance.

  * (Oracle) K-medoid [<cit.>]: number of clusters set to the truth each iteration.
    Data matrix scaled to mean 0 and unit variance.

  * (Oracle) Model based clustering (Mclust) [<cit.>]: number of components set to the truth each iteration.






  * SCAMP: we randomly search for candidate clusters, stopping after we find 50 ·(number of columns). We pre-set α=0.25, m=25, and γ=4 across all simulations. We conduct a single SCAMP iteration per simulation setting and iteration. SCAMP
is provided 16 threads to parallelize the search for candidate clusters.

  * SCAMP20: twenty iterations of the SCAMP procedure are performed and the maximum label heuristic used to cluster the data. Each iteration set to the same parameters as the single SCAMP run.









Observe that methods marked with the prefix “(Oracle)” have had the number of clusters  
provided to them. 
Looking ahead, the subsequent comparison of run-times omits the computation 
cost of estimating the number of clusters.
If the method uses the number of clusters as a tuning parameter, as with k-means,
we provide the number directly.
In the case of affinity propagation, 
parameters are modified to get the method close to the true number.
Model based clustering is capable of making its own estimate of the number 
of components present in a data matrix, 
but must be given a range of possible components to try.
A variety of methods could be used to estimate 
the number of clusters for k-means and k-medoids 
[see, for example, <cit.>, <cit.>, and 
<cit.>].




Results from the simulation study are summarized in figure <ref>.
The first panel displays the distribution of mean ARI over the 50 iterations 
for each of the 32 simulation settings.
There we can see that the methods with an oracle providing the number of clusters 
perform quite well when the simulation parameters 
match their implicit distribution assumptions:
K-means, K-medoids, and Mclust are able to recover cluster assignments almost perfectly in several 
of the multivariate normal settings.

However, there is serious downside risk evident for these same methods: when the parameters of the 
simulation introduce nuisance features and transform observations (through the maps 
(<ref>), described in the appendix),
the cluster assignments can differ substantially from the underlying mixture despite being provided
the true number of mixture components.
On the other hand, we see that SCAMP clusterings usually perform well in terms of ARI,
no matter the simulation settings: it is robust to distributional assumptions.
We also note that the twenty iteration of SCAMP improves the mean ARI under all parameter settings.

For each simulation setting, the method with maximum mean ARI is taken as a baseline value.
The difference between this best baseline ARI and each method's ARI is recorded.
The second panel of figure <ref> shows the distribution of these 
differences across the 32 simulation settings, with a zero value indicating the method performed 
best (according to adjusted rand index) in that setting.
This panel shows that a single SCAMP iteration has ARI within 0.1 of the best performing method 
almost 50% of the time. 
The panel also shows that the twenty SCAMP iterations are almost always performing as well as the 
best performing method, 
and performs the best in terms of adjusted rand index about 25% of the time.

The third panel of figure <ref> shows the distribution of the logarithm 
of mean run-time across methods.
We see that on these simulated data, SCAMP is relatively affordable when given 16 threads: 
a SCAMP single iteration, which estimates the number of clusters, 
is usually finished after 20 seconds.
The run-time comparison may be biased due to implementations: after reviewing the documentation
for the compared methods, we did not see the ability to provide them with additional execution 
threads.



We also see that multiple SCAMP iterations incur additional cost in terms of run-time.

The fourth panel shows the logarithm of the number of clusters estimated by each method 
in the top inset (the methods with an oracle are flat because they are given the truth).
The bottom inset shows the logarithm of the number of clusters with more than 25 observations 
in them.
SCAMP estimate is usually overestimates the truth in both the 10 and 17 cluster scenarios.
Iterating twenty times improves the estimate at the cost of producing an increased number
of smaller clusters.



  §.§.§ Scenario 2: Large sample size, moderate number of clusters


Our second simulation scenario explores settings with larger data matrices: our data matrix increases to 30000 observations.
The graphical overview of figure <ref> 
can be modified to describe this setting by multiplying each cluster size by 10.
We continue to use the 32 settings defined by the parameters described there.



We run a modified collection of methods in this setting due to the computational expense. They are:


  * Leveraged Affinity Propagation (AP)[<cit.>]:
    Using the “apclusterL” implementation in the R package apcluster, 
    we set “maxits=2000”,“convits=2000”, “frac=0.1” (corresponding to 10% of the data), “sweeps=5” and “q=0”.

  * (Oracle) Clustering Large Aplications [<cit.>]: The R implementation clara in the package cluster is used, with k set to the truth, the “samples” parameter set to 50, and sampsize set to 3000, which is 10% of the rows of the data set.

  * (Oracle) K-Means: no change. 

  * (Oracle) K-medoid: no change. 

  * (Oracle) Model based clustering (Mclust):  we still provide the function “Mclust” with the number of clusters. We now initialize EM with 10% of the data randomly selected.



  * SCAMP: no change.

  * SCAMP20: no change.


For each of the 32 simulation settings, we again run 10 iterations per simulation setting.
The results of the simulation are displayed in figure <ref>.
We interpret each panel of the figure as we did in figure <ref>.

In the large matrix setting, we see that SCAMP performs well. 
Iterating SCAMP twenty times continues to improve performance at the cost of additional run-time.
On the other hand, the downside risk of alternative methods is more pronounced in this setting.
For some methods, this deterioration in performance in in part explained by their reliance on 
sub-sampling. Additionally, modifying the parameter settings for Leveraged Affinity Propagation
could potentially improve perfomance in terms of ARI – we kept the current settings
for comparability to scenario 1.





  §.§.§ Scenario 3: Large sample size, large number of clusters


We conclude by simulating data data with a large number of 
observations and a large number of clusters of different sizes.
We continue to generate a data matrix with 30000 rows.
In this scenario, the basic mixture has 
30 components in 20 dimensions, and 56 components in the larger 40 dimensions.
In the 30 cluster composition, cluster sizes range between 100 and 10000 observations.
In the 56 cluster composition, cluster sizes range between 50 and 3950.
Once again, for each of the 32 simulation settings, 
we run 10 iterations.
We refer the reader to the appendix for specific details.

In this scenario, we see that SCAMP continues to perform well. 
Since it does not have to sub-sample,  
it is able to recover much of component structure of the underlying mixture. 
In this scenario, the alternative methods almost all perform worse in terms of ARI than SCAMP.
SCAMP run for twenty iterations almost always performs best in terms of its adjusted rand index. 
A single SCAMP iteration compares well to other methods in terms of run-time,
and slightly over-estimates the true number of clusters on average.







 §.§ Circular Data and Rotated Data


The previous simulation of well-separated mixtures along the measurement axes is the ideal use case for SCAMP: in such cases, it is possible for SCAMP to provide interpretable clusterings of the data matrix. 
Here, we simulate two scenarios where SCAMP will perform poorly. 

Figure <ref> illustrates two examples in ℝ^2. 
In the first case, the data have circular structure. 
SCAMP over partitions these data since the coordinate projections appear highly multimodal, and so the connected structure of the data is lost. 
In the second case, the data are simulated from a mixture of two multivariate Gaussians with mean vectors μ_1 = (1.5,0), μ_2 = (0,1.5) and common covariance matrix 
Σ≡[ 1.0   0.7; 0.7   1.0;           ].
Such data have coordinate projections which appear unimodal. 
As a result, SCAMP clusters the data matrix by producing a single cluster consisting of the entire data matrix. 
The adjusted rand scores reflect SCAMP's difficulty with such data.

If one is willing to sacrifice interpretable clusterings, SCAMP can still be used to cluster 
these kinds of data matrices.
In the first case, SCAMP can be applied to the data matrix taken through 
the map ϕ(x,y) ↦ x^2 + y^2. 
In the second case, SCAMP can be applied to the principal components. 
Figure <ref> shows how SCAMP clusters the data under these transformations. 
Since the clusterings are determined on transformed data matrices, the labels produced by SCAMP 
are not immediately interpretable on the original measurement scale.
However, these clusterings are still useful: SCAMP recovers the underlying number of circular 
components in the first case, the number of mixture components in the second case.
The adjusted rand scores are consequently improved.

In many clustering tasks, a method, such as k-means, is used as a clustering tool after the data 
have been transformed: PCA and kernel methods are two such examples.
These simulations show that in such work flows SCAMP might be an viable alternative. We conclude 
with a practical demonstration of such a work-flow.





 §.§ GTEx Data


As a practical demonstration, we apply SCAMP to bulk RNA seq gene read count data from the 
the Genotype-Tissue Expression (GTEx) Project.
This project is supported by the Common Fund of the Office of the Director of the National Institutes of Health,
and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS.
The data used for the demonstration in this paper were obtained from the
https://www.gtexportal.org/https://www.gtexportal.org/
from the v6p release in the the file “GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_reads.gct.gz”
on April 19, 2018. 
An alternative analysis of the v6 data is given in <cit.>.

This dataset contains per-gene read counts for 56238 genes across 8555 samples.
The samples were collected from 450 human donors.
Each of the 8555 samples is taken from one of 31 primary tissue types from a given donor.
The primary tissue type labels are further refined to 53 tissue labels 
that describe the sample location.
For example 1259 samples are given the primary tissue label “Brain”. 
These 1259 samples have 13 corresponding
tissue labels ranging from “Brain - Amygdala” to  “Brain - Substantia nigra”.
Here, will apply several clustering methods to the unlabeled count data (after transformation).
We will use these two label sets as versions of ground truth for clustering.

These count data are zero inflated and, for a given gene, 
differ by orders of magnitude across the samples.
To apply SCAMP to these data, we first remove 1558 genes 
which have zero counts across all 8555 samples.
This produces a data matrix of integer counts with 8555 
rows and 54680 columns. We transform the counts for
each remaining gene by the map log_2(1+x).
We then normalize each gene across samples by the total sum of the transformed counts.

To apply SCAMP, we next compute the top 50 right-singular 
vectors of the data matrix and then use them rotate the transformed data matrix. 
We then apply SCAMP to the rotated data matrix.
Once again, 
the selection of α can be guided empirically before running SCAMP: 
figure <ref> suggests our default value of α = 0.25 
will label clusters with all SVs exhibiting multimodality (according to the dip test) 
in these data.



For the purpose of demonstration, we will cluster the data with SCAMP 
using all three choices of α since the α values {0.01,0.05,0.25}
correspond to labeling clusters with {12, 13, 14} SVs. 
For each choice of α, we run SCAMP for 100 iterations.
Each search phase of each iteration,
we randomly sample 5000 candidate clusters from the annotation forest.
The final clustering if found using the maximum-vote heuristic across the 100 iterations.

We also cluster the dataset using affinity propagation, Mclust, k-means, and k-medoids.
We modify the parameters for affinity propagation to try to achieve convergence, setting 
“maxits” to 5000, “convits” to 5000, and “lam” to 0.95.
For Mclust, we let the number of mixture components 
range from 2 to 100, and select the final mixture using BIC: in this dataset, 
this selects 24 mixture components.
We use the estimate 24 as the k parameter in k-means and k-medoids.
    


Numerical summaries of the clusterings are listed in the table in figure <ref>.
The different methods perform comparably in terms of VI distance no matter
which set of labels are used as the truth.
k-medoid has the best performance in terms of ARI when the primary 
tissue labels are taken as the truth.  
Affinity propagation has the best performance in terms of ARI when the tissue labels are taken as
the truth.
Visualizing the data can help explain this difference: figure <ref> 
provides a visualization of these data through their t-SNE map.

In figure <ref>, we see both SCAMP 
and affinity propagation both have a tendency to partition primary tissue types into sub-clusters.
As the value of α increases from 0.01 to 0.25, the number of clusters found by SCAMP
increases. This can be seen in the t-SNE visualization, where homogeneous 
islands are split into overlapping clusters across the SCAMP panels.
It would require further analysis to determine if SCAMP and affinity propagation
are detecting groups of biological interest that are not detected by other methods (or each other).
However, since the goal of this section is to show how 
SCAMP can be productively applied to a modern dataset, we defer such analysis to subsequent work.


 
 


§ CONCLUDING DISCUSSION
 

In this paper, we developed a new clustering algorithm called selective clustering annotated using 
modes of projections (SCAMP).
SCAMP relies heavily on the dip test of <cit.>; while developing SCAMP, we also 
developed an extension
of the dip test, described in section <ref>, to test univariate distribution for 
multiple modes.
In section <ref>, we discussed the details of the SCAMP clustering algorithm.
In section <ref>, we showed
that the SCAMP algorithm can used to produce interesting clusterings of many different types of data.

Over the course of the paper, we showed that SCAMP makes minimal distributional assumptions,
has tuning parameters that are relatively easy to set, can produce clusters with interpretable 
labels,
and its able to directly cluster datasets with a large number of observations.
We also discussed one of SCAMP's main limitations:
SCAMP can only detect clusters that are separable along the axes of measurement.
Additionally, we
proposed some work-arounds to this limitation: to apply SCAMP to the GTEx data set, we first
rotated the data through its principal components. 

As our analysis of the GTEx data shows, many clustering algorithms, including SCAMP,
can be used to generate useful clusterings of data matrices with a moderate number of observations.
However, many of these approaches also require the user to provide the number of clusters in the 
data as
a parameter. This requirement necessitates either a separate estimation step, a model comparison 
step, or choosing a large number
for this parameter and then merging clusters. By only requiring that the user set a level α 
for the dip test,
we think SCAMP has simplified the clustering task.

In terms of future research, we are interested in improving SCAMP's annotation step.
Currently SCAMP only annotates clusters relative to the columns in a data set with dip test 
p-values below level α.
This is effectively a variable selection step. 
We think this step causes the SCAMP labels to omit some useful information.
For example, labels do not currently reflect a clusters position relative to a unimodal variables 
with wide shoulders;
it would be useful to improve the annotations to account for such a scenario.

We are also interested the consistency of phase 2 and phase 3 of the SCAMP algorithm.
While there are many hierarchical clustering algorithms, SCAMP's approach of recursively growing 
multiple partition trees and
then selecting leaves across different trees is, to our knowledge, new.
It is of interest to know if there are natural parametric conditions
(or, more generally, separation conditions) for the mixture components 
(<ref>),
under which SCAMP can be shown to recover the number of components in the mixture and to
correctly assign observations to their generating component.




§ ACKNOWLEDGEMENTS


EG wishes to thank both K. K. Dey and V. Voillet for their separate discussions of the GTEx data. 
This work was funded by a grant from the Bill and Melinda Gates Foundation to RG [OPP1032317]
and a grant from the NIGMS to GF [R01 GM118417-01A1].









§ APPENDIX
 

In these appendices, we provide proofs for section <ref> 
and details for the separated mixture simulation in section <ref>.



 §.§ Proofs
 

In this appendix, we provide proofs for lemmas and theorems found in section <ref>.


Proof of Lemma <ref>


   Let q_α denote the α-quantile of the F-distribution, so that F(q_α-) ≤α and
   F(q_α) ≥α. For a fixed value of N, select the 2N-1 quantiles
  
    q_1≡ q_1/2N, q_2 ≡ q_2/2N,…, q_2N-1≡ q_2N-1/2N .

  Suppose first that F is continuous. This means the values (<ref>) are distinct.
  Begin by considering H_1 ≡ F/F(q_1-). Per the discussion starting section 5 of <cit.>.,
  there exists a convex distribution function C_1 on (-∞,q_1) such that
  
    sup_x<q_1 |H_1-C_1| ≤ 1/2  .

  Define G_1 ≡ F(q_1-)· C_1. This is convex on (-∞,q1). Then
  
    sup_x < q_1 |F-G_1| = F(q_1-) ·sup_x < q_1 |F/F(q_1-)-C_1|≤1/4N .

  Next, define H_N ≡ [F(x)-F(q_2N-1)]/[1-F(q_2N-1)]. This gives a distribution function on [q_2N-1,∞).
  There exists a concave distribution function C_N on [q_2N-1,∞) such that
  
    sup_x ≥ q_2N-1 |H_N-C_N| ≤ 1/2  .

  Define G_N ≡ F(q_2N-1)+(1-F(q_2N-1))C_N, which is concave. Then
  
    sup_x ≥ q_2N-1 |F-G_N| = (1-F(q_2N-1)) ·sup_x > q_2N-1 |H_N-C_N|≤1/4N .

  Finally, for j ∈{2,3,…,2N-1}, define
  
    H_j ≡[F(x)-F(q_j-1)]/F(q_j-)-F(q_j-1)

  on [q_j-1,q_j). The uniform d.f. on the interval [q_j-1,q_j], which we denote by C_j, has the trivial bound
  
    sup_q_j-1≤ x  ≤ q_j |H_j-C_j| ≤ 1  .

  Define G_j ≡ F(q_j-1)+(F(q_j-)-F(q_j-1))C_j. Then
  
    sup_q_j-1≤ x < q_j |F-G_j| = (F(q_j-)-F(q_j-1)) sup_q_j-1≤ x < q_j |H_j-C_j|≤ F(q_j-)-F(q_j-1) ≤1/2N .

  Finally define
  
    G(x) ≡∑_j=1^2N G_j 1_[q_j-1≤ x < q_j) ,

  where we define q_0 ≡ -∞ and q_2N = ∞. The linear components of G(x) are both convex and concave,
  so by construction G(x) ∈𝒰_N and ρ(F,G) ≤ 1/2N.

  If F is not continuous, the quantiles (<ref>) may no longer be distinct.
  However, the construction used in the continuous case still essentially works. Define
  
    L ≡sup{x | F(x) ≤1/2N}

  and
  
    U ≡inf{x | F(x) ≥2N-1/2N} .

  Then the constructions at (<ref>) and (<ref>) still produce the approximations on (-∞,L)
  and [U,∞) that are smaller than 1/4N.

  If L=U, then to conform with 𝒰_N we may proceed as follows: pick ϵ  > 0 and
  apply the construction (<ref>) on the interval [U+ϵ,∞). Then, divide the interval [U,U+ϵ)
  equally according to the choice of N, and apply the linear approximation of (<ref>) to each component. Subsequently
  define G(x) in (<ref>) using these approximations.

  If L < U, then there are distinct quantile values (<ref>). Repeated value in the list corespond to jump points,
  and we may repeat the method of the case L=U:
  insert linear approximations after jump points as needed to conform with 𝒰_N.


Proof of Lemma <ref>


  When N = 1, this is proved by Theorem 1 of <cit.>. For N > 1, the proof 
  of Hartigan and Hartigan extends naturally. Suppose N > 1.
  For G ∈𝒰_N, define
  
    H(x) ≡ G(0)1_[x < 0]+G(x)1_[0 ≤ x ≤ 1] +G(1)1_[x > 1] .

  Recalling definition <ref>, suppose first m_N < 0. Then setting
  
    m_H_1 =c_H_1 = m_H_2 = … = m_H_N-1 = c_H_N-1 =  m_H_N≡ 0

  places H(x) ∈𝒱_N. If m_1 > 1, setting all modal and antimodal values to 1 places H(x) ∈𝒱_N
  Similarly, in the intermediate cases where m_j < 0, c_j > 0 , c_k < 1, and m_k+1 > 1 with j < k, setting
  all modal and antimoal values below 0 to 0 and all values above 1 to 1 places H(x) ∈𝒱_N. Thus,
  
    ρ(F,H) =  sup_x |F(x)-H(x)| = sup_0 ≤ x ≤ 1 |F(x)-G(x)| ≤ρ(F,G) .

  Hence ρ(F,𝒱_N) ≤ρ(F,𝒰_N).

  Conversely, pick G ∈𝒱_N. By Lemma <ref>, it suffices to consider only G such that
  ρ(F,G) ≤ 1/2N. Such G satify 0 ≤ G(0) ≤ 1/2N < 1-1/2N ≤ G(1) ≤ 1.
  Define α≡ 0 ∨ G(0) , β≡ 1 ∧ G(1), and
  
    H(x) ≡α 1_[G < α] +G(x) 1_[α≤ G(x) ≤β] + β 1_[G > β] .

  Since G is non-decreasing on [0,1] by definition, so too is H.
  Define L ≡sup{x | G(x) ≤α} and U ≡inf{x | G(x) ≥β}.
  For the values m_1,c_1,…,c_N-1,m_N that fall below α, set all of them to G(L-) if the largest is a point of
  convexity. If the largest is a point of concavity, set that value to G(L) (capturing the jump) and the remainder to G(L-).
  Perform the analogous procedure for those points above β. Doing so, we see H(x) ∈𝒱_N, and
  
    ρ(F,H) = sup_α≤ G(x) ≤β|F(x)-G(x)| ≤ρ(F,G) .

  Now for a ≥ 1 define
  
    G_a(x) ≡
          GCM(H· 1_[x ≥ -a], (-∞,m_1)) 
     
          H(x)                         for  x ∈ [m_1,m_N] 
    
          LCM(H· 1_[x ≤ a], (m_N,∞)) 
    
         .

  Here G_a(x) ∈𝒰_N.  This extends the G_a of Theorem 1 in <cit.>. The remainder
  of the proof follows the argument in Hartigan and Hartigan, with the expected notational changes. We include it for reference.

  
  The function G_a is constant on (-∞,a], consists of a line segment on [-a,x_1] for
  some 0 ≤ x_1 ≤ m_1, and is equal to H on [x_1,m_1]. Here we also have that as a ↗∞,
  the slope of the linear segment tends to zero. Since H is non-decreasing,
  
    sup_0 ≤ x ≤ m_1 |H(x)-G_a(x)| → 0  as  a ↗∞ .

  By symmetry
  
    sup_m_N ≤ x ≤ 1 |H(x)-G_a(x)| → 0  as  a ↗∞ .

  Since G_a(x) = H(x) for x ∈ [m_1,m_N], we conclude as a ↗∞ that
  
    ρ(F,G_a) = sup_0 ≤ x ≤ 1  |G_a(x) - F(x)| →sup_0 ≤ x ≤ 1|H(x) - F(x)| = ρ(F,H) .

  Thus for each ϵ > 0 and G∈𝒱_N, there exists G_a ∈𝒰_N with
  
    ρ(F,G_a) ≤ρ(F,G) + ϵ .

  Thus ρ(F,𝒰_N) ≤ρ(F,𝒱_N).


Proof of Lemma <ref>


  The argument of Theorem 2 in <cit.> applies here, with D_N replacing D, and
  Lemma <ref> of this paper replacing Theorem 1 in <cit.>.


Proof of Theorem <ref>


  The argument of Theorem 3 in <cit.> applies here, with D_N replacing D, and
  Lemma <ref> of this paper replacing Theorem 2 in <cit.>.




Proof of Theorem <ref>


  As noted in Section <ref>, the claim essentially follows from Theorem 5 of <cit.>.
  We recapitulate the argument of Hartigan and Hartigan here, providing 
  details showing that the conditions (<ref>) through
  (<ref>) imply
  
    sup_-∞≤ x_1 ≤ x_2 ≤∞[F(x_2)-F(x_1)]^k/|F(x_1)+F(x_2) - 2F(x)| < ∞

  so long as
  
    m_1,…,m_N,c_1,…,c_N-1∉ (x_1,x_2) .

  Once demonstrated, the argument given in the second half of Theorem 5 in <cit.> proves the claim. We need to consider new details, since the additional peaks and valleys of multimodality mean we must consider
  some additional, but unsurprising, cases. 

  Without loss of generality, suppose m_1 = 0. To begin, if x_1,x_2 < m_1=0 or x_1,x_2 > m_N
  (and so (<ref>) or (<ref>) hold), (<ref>) is proved by the argument given in the first half of Theorem 5 in <cit.>
  We now consider condition (<ref>). In the following analysis of cases, fix j between 1 ≤ j ≤ N.

  Case 1
  
  Pick c_j-1+ϵ_1 < x_1, x_2 < m_j -ϵ_2,  F'(x) > 0. Then,
  
    d/dxlog F'(x) > B ≡ B_ϵ_1∧ B_ϵ_2 > 0

  Thus
  
    ∫_x_1^x_2d/dxlog F'(x) = logF'(x_2)/F'(x_1) > B(x_2-x_1) ,

  and so
  
    F'(x_2)/F'(x_1)≥exp[B(x_2-x_1)]  .

  This implies
  
    F(x_2)-F(x_1+x_2/2) ≥[F(x_1+x_2/2)-F(x_1)]exp[B(x_2-x_1/2)] .

  To see that this is so, for y ∈ [(x_1+x_2)/2,x_2] ≡ℐ, define the functions
  
    h(y)    ≡ F(y) , 
    
        g(y)    ≡ F(y-x_1+x_2/2+x_1)  .

  
  Then by Cauchy's mean value theorem, there exists a c ∈ ([x_1+x_2]/2,x_2) such that
  
    h(x_2) - h(x_1+x_2/2)/g(x_2) - g(x_1+x_2/2) = h'(c)/g'(c) .

  By their definition, this is equivalent to
  
    F(x_2) - F(x_1+x_2/2)/F(x_1+x_2/2)-F(x_1) = F'(c)/F'(c-x_1+x_2/2+x_1)≥exp(B[x_1+x_2/2-x_1])
        =exp(B[x_2-x_1/2]) ,

  with the inequality following by (<ref>), justifying the claim. Let x≡ (x_1+x_2)/2 and δ≡ (x_2-x_1)/2. Observing Bδ≥ 0, we see (<ref>) implies
  
    1/1-[F(x)-F(x_1)]  / [F(x_2)-F(x)]≤1/1-exp(-B δ)

  and
  
    [1 + F(x)-F(x_1)/F(x_2)-F(x)]^k ≤[1+exp(-Bδ)]^k ≤ 2^k .

  Now if B δ≥ 1,  we combine the preceding to find
  
    [F(x_2)-F(x_1)]^k/|F(x_2)+F(x_1) - 2 F(x)|≤[F(x_2)-F(x)]^k-1 2^k/1-exp(-B δ)≤2^k/1-1/e < ∞ .

  Since 2/x > 1/[1-exp(-x)] for 0 < x < 1, if 0 < Bδ < 1 we find
  
    [F(x_2)-F(x_1)]^k/|F(x_2)+F(x_1) - 2 F(x)|≤[F(x_2)-F(x)]^k-1 2^k/1-exp(-B δ)≤[F(x_2)-F(x)]^k-12^k+1/Bδ .

  The mean value theorem gives for some η∈ ([x_1+x_2]/2,x_2) that
  
    F(x_2)-F(x)=F(x_2)-F(x_1+x_2/2) = x_2-x_1/2F'(η) = δ F'(η) .

  At the mode m_j, we also have F'(η) ≤ F'(m_j). We thus continue from (<ref>) to see
  
    [F(x_2)-F(x_1)]^k/|F(x_2)+F(x_1) - 2 F(x)|≤δ^k-2[F'(m_j)]^k-12^k+1/B  < [F'(m_j)]^k-12^k+1/B^k-1  < ∞ .


  Case 2
  
  Since we assume there is a k ≥ 2 for which F^(k)(m_j) ≠ 0, pick ϵ small
  enough that the sign of F^(k)(x) does not change for x ∈ [m_j-ϵ,m_j]
  and so that F^(x)(x) ≠ 0. 
  Now suppose m_j-ϵ < x_1 < x_2 < m_j.   
  Two applications of Taylor's theorem give
  
  
  
  
  
  
  
  
    F(x_2)+F(x_1)-2F(x) = (x_2-x_1)^2/8[F”(y_1)+F”(y_2)] .

  with  y_1 ∈ (x_1,x) and y_2 ∈ (x,x_2). 
  If, excepting the first derivative, F^(k)(m_j) is the first non-zero derivative at m_j, the Taylor polynomial of F”(y_1) is
  
    P_k-3(y_1-m_j) = F”(m_j)+F”'(m_j)(y_1-m_j)+…+ F^(k-1)(m_j)/(k-3)! (y_1-m_j)^k-3 = 0 ,

  with remainder 
  
    R_k-3(y_1-m_j) = F^(k)(η_y_1)/(k-2)! (y_1-m_j)^k-2 .
      
  The analogous expansion of F”(y_2) at m_j gives
  
    |F(x_2)+F(x_1)-2F(x)|    =
      (x_2-x_1)^2/8(k-2)!· |F^(k)(η_y_1)(y_1-m_j)^k-2 + F^(k)(η_y_2)(y_2-m_j)^k-2|
      
       ≥(x_2-x_1)^2/8(k-2)!· |F^(k)(η_k)(y_1-m_j)^k-2 + F^(k)(η_k)(y_2-m_j)^k-2|
      
       ≥(x_2-x_1)^2 · |F^(k)(η_k)|· |(m_j-x)^k-2|/8(k-2)!

  with η_k ∈{η_y_1,η_y_2} chosen as a function of k to minimize
  the sum (<ref>).
  
  
  
  
  Again using the mean value theorem, we have
  
    F(x_2)-F(x_1) = (x_2-x_1)F'(ζ) ≤ (x_2-x_1)F'(m_j) .

  Since (x_2-x_1)/(m_j-x) ≤ 2, we find
  
    [F(x_2)-F(x_1)]^k/|F(x_2)+F(x_1)-2F(x)|≤8(x_2-x_1)^k-2[F'(m_j)]^k · (k-2)!/|F^(k)(η_k)|· |m_j-x|^k-2≤ 2^k+1[F'(m_j)]^k· (k-2)!/|F^(k)(η_k)| < ∞ .


  Case 3

  Next suppose c_j-1 < x_1 < x_2 < c_j-1+ϵ. A similar analysis to that given in case 2 shows
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    |F(x_2)+F(x_1)-2F(x)|
        
        ≥(x_2-x_1)^2|F^(k)(η_k)|· |(x-c_j-1)^k-2|/8(k-2)! .

  
  
  
  
  However, the mean value theorem and increase in x still gives
  
    F(x_2)-F(x_1) = (x_2-x_1)F'(ζ) ≤ (x_2-x_1)F'(m_j) ,

  and so the analogous bound holds.
  
  
  
  

  Case 4

  Now suppose c_j-1+ϵ_1 < x_1 ≤ m_j-ϵ_2 < m_j-η≤ x_2 < m_j, and that F(m_j)-F(m_j-η) ≤ F(m_j-η) - F(m_j-ϵ_2). The assumptions imply:
  
    F(x_2)-F(m_j-η) ≤ F(m_j)-F(m_j-η) ≤ F(m_j-η)-F(m_j-ϵ_2) ≤ F(m_j-η)-F(x_1) .

  From this it follows
  
    F(m_j-η) + [F(m_j-η) - F(x_1) ] - F(x_2) ≥ F(m_j-η) + [F(m_j) - F(m_j-η)] - F(x_2) = F(m_j)-F(x_2) ≥ 0 .

  Re-arranging, this implies
  
    F(x_2)-F(x_1) ≤ 2[F(m_j-η)-F(x_1)] .

  Re-arranging and using convexity, it also implies
  
    F(m_j-η) ≥F(x_1)+F(x_2)/2≥ F(x) .

  
  
  
  
  
  Now define
  
    h_x_1(t) ≡ F(t) + F(x_1) - 2F(t+x_1/2) ,

  with t ∈ (x_1, m_j). Recall F is convex and non-decreasing on [c_j-1+ϵ_1,m_j]. Differentiating, we thus see
  
    h_x_1'(t) = F'(t)- F'(t+x_1/2)  > 0 .

  Since m_j-η < x_2, we thus have
  
    h_x_1(m_j-η) ≤ h_x_1(x_2) ⟺ F(m_j-η) + F(x_1) - 2F(m_j-η+x_1/2) ≤ F(x_2) + F(x_1) - 2F(x_2+x_1/2) .

  This is equivalent to
  
    1/F(x_2) + F(x_1) - 2F(x_2+x_1/2)≤1/F(η) + F(x_1) - 2F(m_j-η+x_1/2) .

  Combining this with (<ref>), we infer
  
    [F(x_2)-F(x_1)]^k/F(x_2) + F(x_1) - 2F(x_2+x_1/2)≤2^k[F(m_j-η)-F(x_1)]^k/F(m_j-η) + F(x_1) - 2F(m_j-η+x_1/2)≤ 2^kC_η < ∞ ,

  with the penultimate inequality following by Case 1.

  Case 5
  
   Suppose now that c_j-1 < x_1 ≤ c_j-1+η < c_j-1+ϵ_1 ≤ x_2 ≤ m_j-ϵ_2, and that F(c_j-1+η)-F(c_j-1) ≤ F(c_j+1+ϵ_1) - F(c_j-1+η). Similar to Case 4,
  
    F(c_j-1+η)-F(x_1) ≤ F(c_j-1+η)-F(c_j-1) ≤ F(c_j-1+ϵ_1)-F(c_j-1+η) ≤ F(x_2)-F(c_j-1+η) .

  From this follows
  
    F(x_2)-F(x_1) ≤ 2[F(x_2)-F(c_j+1+η)] .

  A symmetric argument with Case 4 gives
  
    [F(x_2)-F(x_1)]^k/F(x_2) + F(x_1) - 2F(x_2+x_1/2)≤2^k[F(x_2)-F(c_j-1+η)]^k/F(x_2) + F(c_j-1+η) - 2F(c_j-1+η+x_2/2)≤ 2^k C_η < ∞ ,


  Case 6
  
  Finally, suppose
  
    c_j-1 < x_1 ≤ c_j-1+η_1 ≤  c_j-1+ϵ_1  < m_j-ϵ_2 ≤ m_j - η_2 ≤ x_2 < m_j  ,

  with both
  
    F(c_j-1+η_1)-F(c_j-1) ≤ F(c_j+1+ϵ_1) - F(c_j-1+η_1)

  and
  
    F(m_j)-F(m_j-η_2) ≤ F(m_j-η_2) - F(m_j-ϵ_2) .

  Then
  
    F(c_j-1+η_1)-F(x_1) ≤ F(c_j-1+η_1)-F(c_j-1) ≤ F(c_j+1+ϵ_1) - F(c_j-1+η_1) ≤ F(x_2) - F(c_j-1+η_1)

  and
  
    F(x_2)-F(m_j-η_2) ≤ F(m_j)-F(m_j-η_2) ≤ F(m_j-η_2)-F(m_j-ϵ_2) ≤ F(m_j-η_2)-F(x_1) .

  Jointly these imply
  
    F(x_2)-F(x_1) ≤ 2[F(m_j-η_2)-F(c_j-1+η_1)] .

  Combining the argument of cases 4 and 5 gives the bound
  
    [F(x_2)-F(x_1)]^k/F(x_2) + F(x_1) - 2F(x_2+x_1/2)≤2^k[F(m_j-η_2)-F(c_j-1+η_1)]^k/F(m_j-η_2) + F(c_j-1+η_1) - 2F(c_j-1+m_j+η_1-η_2/2)≤ 2^k C_η_1,η_2 < ∞ .


  Conclusion

  Taken together, cases 1-6 show condition (<ref>) implies (<ref>). Symmetric arguments imply when (<ref>) holds,
  (<ref>) holds. We thus see (<ref>) follows from the given assumptions.
  
  To show √(n)D_N(F_n) →_p 0, construct G_n as follows: let G_n be the LCM of F_n in [-∞,m_1], the GCM of F_n on [m_N,∞], the LCM of
  F_n on intervals [c_j-1,m_j], and the GCM of F_n on intervals [m_j,c_j]. Then G_n ∈𝒰_N. 
  Considering maximal linear components of G_n that occurs between the atoms m_1,…,m_N,c_1,…,c_N-1, such that G_n=F_n at the end-points. The argument of Hartigan and Hartigan in Theorem 5 of <cit.>
  show these intervals are of length o_P(1). Convergence in probability follows.






 §.§ Separated Mixture Settings
 



This section of the appendix describes our settings for the separated mixtures simulation.

When all parameters are set to their default, the mixture components of (<ref>) are all multivariate Gaussian.
The mean vector of each component is randomly selected so its entries lie in {0,6}.
The variance-covariance matrix is randomly selected from the space of positive definite matrices with variance bounded above by 3.


We now describe how parameter settings affect the simulation.
When the “Second Mixture Dimension” parameter is set to 20, a random sample from a 
second mixture of either size 3000 or 30000 is taken. The number of rows 
depends on the simulation scenario.




This leads to a data matrix of size 3000 × 40 or 30000 × 40, depending on the 
simulation scenario.
The clustering ground truth is taken to be the cross-product of the labels of the mixture
components which produces a row of the data matrix.
This creates 17 distinct clusters ranging in size between 125 and 375 observations.

The “Second Mixture Mean Different” parameter only affects the simulation if the “Second Mixture
Dimension” parameter is set to 20.
If the “Second Mixture Mean Different” parameter is set to “False”, the mean 
vectors of components in the simulated mixture have entries randomly chosen from
{0,6}.
If the “Second Mixture Mean Different” parameter is set to “True”, the
mean vectors of components in the simulated mixture have entries randomly chosen from
{0,3,6}.



If the “Noise Component Dimension” parameter is set to 20, a sample of 3000 
or 30000 observations is taken from a multivariate t distribution with 5 degrees of freedom.
The non-centrality parameter is set to 2.5 for each component.
The scale matrix sampled from the space of positive definite matrices with diagonal elements 
bounded above by 3.
This sample is adjoined to the samples from the mixture distribution(s). It is mean to represent 
the case when measurements are taken on variables unrelated to the signal of interest.
Since the sample is from a single distribution, it does not affect the clustering ground truth.

If the “Transform Coordinates” parameter is set to “True”, rows sampled from the mixtures are taken through the following maps:

    f_1(x)    = (√(|x_1|),…,√(|x_p|)) . 
    
      f_2(x)    = (exp(x_1),…,exp(x_p)) . 
    
      f_3(x)    = (x_1^2,…,x_p^2) . 
    
      f_4(x)    = (x_1,…,x_p) .

The maps are sequentially applied to the components of the mixture, starting with f_1.
Once all four maps have been applied to different mixture components, the process repeats
from f_1, until all clusters have been taken through a map.

Finally, if the “Sampling Regime” parameter is set to “T plus Gaussian”, components of the mixture alternate between multivariate normal and multivariate t with 5 degrees of freedom.
For the multivariate t, the non-centrality parameter and scale matrix are sampled from the same space as its Gaussian counterpart.
If the “Sampling Regime” parameter is set to “Gaussian”, components of the mixture are all initially samples from multivariate normal.



  §.§.§ Scenario 1: Moderate sample size, moderate number of mixture components


To simulate this scenario, we begin by taking a sample of size 3000 from a mixture in 20 dimensions.
The mixture components have the following sizes, and are sampled in the stated order:

    Mixture Component Sizes ≡ {1000,500,250,250,250,250,125,125,125,125} .

In this setting when the “Second Mixture Dimension” parameter is set to 20, a random sample from a second mixture of size 3000 is taken.
The mixture components have the following sizes; samples are taken the stated order:

    Second Mixture Component Sizes ≡ {250,125,250,125,125,250,500,125,250,1000} .

As before, this creates 17 distinct clusters.



  §.§.§ Scenario 2: Large sample size, moderate number of mixture components


The mixture components have the following sizes, and are sampled in the stated order:

    Mixture Component Sizes ≡ {10000,5000,2500,2500,2500,2500,1250,1250,1250,1250} .

In this setting when the “Second Mixture Dimension” parameter is set to 20, a random sample from a second mixture of size 30000 is taken.
The mixture components have the following sizes; samples are taken the stated order:

    Second Mixture Component Sizes ≡ {2500,1250,2500,1250,1250,2500,5000,1250,2500,10000} .

As before, this creates 17 distinct clusters.



  §.§.§ Scenario 3: Large sample size, moderate number of mixture components


The mixture components have the following sizes, and are sampled in the stated order:

    Large Sample Mixture Component Sizes ≡{   10000,5000,2500,2500,1250,1250
       1250,1250,500,500,500,500
       250,250,250,250,250,250
       250,250,100,100,100,100
       100,100,100,100,100,100
    } .

In this setting when the “Second Mixture Dimension” parameter is set to 20, a random sample from a second mixture of size 30000 is taken.
The mixture components have the following sizes; samples are taken the stated order:

    Large Sample Second Mixture Component Sizes ≡{   250,100,250,100,100,250
       500,100,100,100,1250,250
       500,100,250,1250,100,500
       5000,250,10000,2500,250,1250
       100,100,500,2500,1250,250
    } .

When the mixtures (<ref>) and (<ref>) are combined, they 
produce a mixture with 56 components, with sample sizes
oranging between 50 and 3950 observations.



authordate1


