Modal clustering asymptotics with applications to bandwidth selection

TL;DR

This paper investigates the asymptotic properties of density estimation in modal clustering, focusing on bandwidth selection to accurately identify clusters based on density modes.

Contribution

It introduces an asymptotic framework for density estimation in modal clustering and proposes a method for bandwidth selection based on this theory.

Findings

Asymptotic approximation of the partition distance metric

Guidelines for bandwidth selection in modal clustering

Improved understanding of density estimation impact on clustering accuracy

Abstract

Density-based clustering relies on the idea of linking groups to some specific features of the probability distribution underlying the data. The reference to a true, yet unknown, population structure allows to frame the clustering problem in a standard inferential setting, where the concept of ideal population clustering is defined as the partition induced by the true density function. The nonparametric formulation of this approach, known as modal clustering, draws a correspondence between the groups and the domains of attraction of the density modes. Operationally, a nonparametric density estimate is required and a proper selection of the amount of smoothing, governing the shape of the density and hence possibly the modal structure, is crucial to identify the final partition. In this work, we address the issue of density estimation for modal clustering from an asymptotic perspective. A natural and easy to interpret metric to measure the distance between density-based partitions is discussed, its asymptotic approximation explored, and employed to study the problem of bandwidth selection for nonparametric modal clustering.

Figures

Tables

Table 1. Table 1: Top panel: the EDM (solid line), the AEDM (dashed grey line), and the bounds AB1 (dotted line) and AB2 (dot-dashed line) versus h ℎ h , for n = 100 , 1000 , 10000 𝑛 100 1000 10000 n=100,1000,10000 . All the expressions are evaluated by assuming f 𝑓 f and all the involved quantities known. The minimum EDM is reported below the plots, together with the EDM for the oracle bandwidths h AEDM subscript ℎ AEDM h_{\rm AEDM} and h MISE , 1 subscript ℎ MISE 1 h_{\rm MISE,1} . Bottom panel: average distances in measure (and their standard error) for the proposed bandwidth selectors and the plug-in bandwidth for density gradient estimation. Results refer to density M1.

	n =100	n=1000	n=10000
$h_{\mathrm{EDM}}$	0.144	0.060	0.020
$h_{\mathrm{AEDM}}$	0.164	0.103	0.050
$h_{\mathrm{MISE},1}$	0.146	0.081	0.044
${\hat{h}}_{\mathrm{AEDM}}$	0.331 (0.148)	0.150 (0.164)	0.034 (0.076)
${\hat{h}}_{\mathrm{AB1}}$	0.297 (0.161)	0.111 (0.134)	0.027 (0.060)
${\hat{h}}_{\mathrm{AB2}}$	0.320 (0.150)	0.128 (0.147)	0.028 (0.060)
${\hat{h}}_{\mathrm{PI},1}$	0.221 (0.176)	0.063 (0.084)	0.029 (0.052)

Table 2. Table 2: Cf. Table 1 . Results refer to density M2.

	n =100	n=1000	n=10000
$h_{\mathrm{EDM}}$	0.131	0.040	0.008
$h_{\mathrm{AEDM}}$	0.143	0.047	0.008
$h_{\mathrm{MISE},1}$	0.165	0.041	0.011
${\hat{h}}_{\mathrm{AEDM}}$	0.437 (0.195)	0.201 (0.122)	0.009 (0.005)
${\hat{h}}_{\mathrm{AB1}}$	0.382 (0.201)	0.185 (0.118)	0.008 (0.005)
${\hat{h}}_{\mathrm{AB2}}$	0.423 (0.197)	0.196 (0.118)	0.008 (0.005)
${\hat{h}}_{\mathrm{PI},1}$	0.256 (0.159)	0.092 (0.076)	0.008 (0.005)

Table 3. Table 3: Cf. Table 1 . Results refer to density M3.

	n =100	n=1000	n=10000
$h_{\mathrm{EDM}}$	0.045	0.010	0.003
$h_{\mathrm{AEDM}}$	0.051	0.016	0.011
$h_{\mathrm{MISE},1}$	0.054	0.034	0.022
${\hat{h}}_{\mathrm{AEDM}}$	0.071 (0.072)	0.042 (0.086)	0.017 (0.036)
${\hat{h}}_{\mathrm{AB1}}$	0.059 (0.059)	0.038 (0.075)	0.017 (0.031)
${\hat{h}}_{\mathrm{AB2}}$	0.067 (0.067)	0.040 (0.080)	0.017 (0.033)
${\hat{h}}_{\mathrm{PI},1}$	0.050 (0.072)	0.024 (0.025)	0.019 (0.017)

Table 4. Table 4: Cf. Table 1 . Results refer to density M4.

	n =100	n=1000	n=10000
$h_{\mathrm{EDM}}$	0.039	0.009	0.003
$h_{\mathrm{AEDM}}$	0.187	0.040	0.005
$h_{\mathrm{MISE},1}$	0.040	0.015	0.005
${\hat{h}}_{\mathrm{AEDM}}$	0.066 (0.088)	0.018 (0.054)	0.011 (0.037)
${\hat{h}}_{\mathrm{AB1}}$	0.057 (0.077)	0.014 (0.032)	0.009 (0.033)
${\hat{h}}_{\mathrm{AB2}}$	0.063 (0.082)	0.015 (0.035)	0.009 (0.033)
${\hat{h}}_{\mathrm{PI},1}$	0.051 (0.069)	0.011 (0.014)	0.005 (0.005)

Table 5. Table 5: Cf. Table 1 . Results refer to density M5.

	n =100	n=1000	n=10000
$h_{\mathrm{EDM}}$	0.058	0.012	0.005
$h_{\mathrm{AEDM}}$	0.059	0.012	0.005
$h_{\mathrm{MISE},1}$	0.058	0.013	0.005
${\hat{h}}_{\mathrm{AEDM}}$	0.222 (0.184)	0.015 (0.018)	0.005 (0.003)
${\hat{h}}_{\mathrm{AB1}}$	0.195 (0.180)	0.013 (0.010)	0.005 (0.003)
${\hat{h}}_{\mathrm{AB2}}$	0.212 (0.184)	0.014 (0.011)	0.005 (0.003)
${\hat{h}}_{\mathrm{PI},1}$	0.179 (0.158)	0.013 (0.009)	0.005 (0.003)

Table 6. Table 6: Minimum EDM associated with a density estimate with bandwidth matrix 𝐇 𝐇 \mathbf{H} selected to minimize the EDM ( 𝐇 EDM subscript 𝐇 EDM \mathbf{H}_{\rm EDM} ) and the MISE for gradient estimation ( 𝐇 MISE , 1 subscript 𝐇 MISE 1 \mathbf{H}_{\rm MISE,1} ). Different parametrizations for 𝐇 𝐇 \mathbf{H} are considered. In both cases, the true density as well as all the involved quantities are assumed to be known. Results refer to density M6.

	${𝐇}_{\mathrm{EDM}}$		${𝐇}_{\mathrm{MISE},1}$
	n =100	n=1000	n =100	n=1000
	0.006	0.004	0.064	0.040
	0.006	0.004	0.064	0.040
	0.005	0.003	0.042	0.024

Table 7. Table 7: Cf. Table 6 . Results refer to density M7.

	${𝐇}_{\mathrm{EDM}}$		${𝐇}_{\mathrm{MISE},1}$
	n =100	n=1000	n =100	n=1000
	0.114	0.044	0.116	0.054
	0.114	0.042	0.115	0.055
	0.110	0.040	0.121	0.054

	Components	${\pi }_k$	${\mu }_k$	${\sigma }_k^2$
	1	0.75	0.00	0.83
	2	0.25	1.37	0.09

	Components	${\pi }_k$	${\mu }_k$	${\Sigma }_k$
	1	0.5	$\left(\begin{array}{c}1\\ -1\end{array}\right)$
	2	0.5	$\left(\begin{array}{c}-1\\ 1\end{array}\right)$

	Components	${\pi }_k$	${\mu }_k$	${\Sigma }_k$
	1	0.43	$\left(\begin{array}{c}-1\\ 0\end{array}\right)$
	2	0.43	$\left(\begin{array}{c}1\\ 1.15\end{array}\right)$
	3	0.14	$\left(\begin{array}{c}1\\ -1.15\end{array}\right)$

Equations

$${\nu }_x^{\prime }(t)=-\nabla f({\nu }_x(t)),{\nu }_x(0)=x.$$

$$W_{-}^u(x_0)=\{x\in {\mathbb{R}}^d:\underset{t\to -\infty }{\lim }{\nu }_x(t)=x_0\}.$$

$$x^{(l+1)}=x^{(l)}+A\frac{\nabla f(x^{(l)})}{f(x^{(l)})},$$

$${\hat{f}}_h(x)=\frac{1}{nh}\sum _{i=1}^{n}K\left(\frac{x-X_i}{h}\right),$$

$$\mathrm{ISE}(h)={\int }_{\mathbb{R}}\{{\hat{f}}_h(x)-f(x){\}}^{2}dx.$$

$$\mathrm{MISE}(h)=\mathbb{E}\left[\mathrm{ISE}(h)\right]$$

$$d(\mathcal{C},\mathcal{D})=\frac{1}{2}\underset{\sigma \in {\mathcal{P}}_s}{\min }\left\{\sum _{i=1}^r\mathbb{P}(C_i\Delta D_{\sigma (i)})+\sum _{i=r+1}^s\mathbb{P}(D_{\sigma (i)})\right\},$$

$${\rm EDM}(h)=\mathbb{E}\big{[}d(\hat{\mathscr{C}}_{h},\mathscr{C}_{0})\big{]},$$

$$\psi (\mu ,{\sigma }^2)$$

$${\rm AEDM}(h)=\sum_{j=1}^{r-1}\frac{f(m_{j})}{f^{(2)}(m_{j})}\psi\Big{(}\tfrac{1}{2}\mu_{2}(K)f^{(3)}(m_{j})h^{2},R(K^{(1)})f(m_{j})(nh^{3})^{-1}\Big{)},$$

$$\psi (\mu ,{\sigma }^2)\le (2/\pi {)}^{1/2}\sigma +|\mu |$$

$$\mathrm{AB}1(h)$$

$$\mathrm{AB}2(h)$$

$$+(32\pi {)}^{-1/2}{\mu }_2(K{)}^{2}R(K^{(1)}{)}^{-1/2}{a}_2{n}^{1/2}{h}^{11/2}.$$

$$a_{j1}$$

$$a_{j2}$$

$$h_{AB1}$$

$$h_{AB2}$$

$$h_{\mathrm{AMISE},1}={\left(\frac{3R(K^{(1)})}{{\mu }_2(K{)}^{2}R(f^{(3)})}\right)}^{1/7}{n}^{-1/7},$$

$${\hat{h}}_{\mathrm{crit}}=\inf \{h>0:{\hat{f}}_h(\cdot )\text{mbox}hasexactlyonemode\},$$

$${\hat{f}}_{\mathbf{H}}(\mathbf{x})=\frac{1}{n}\sum _{i=1}^n|\mathbf{H}{|}^{-1/2}K\left({\mathbf{H}}^{-1/2}(\mathbf{x}-{\mathbf{X}}_i)\right),$$

$$\mathrm{EDM}(h)=\sum _{j=1}^{r-1}\mathbb{E}|F({\hat{m}}_{h,j})-F(m_j)|.$$

$$F(\widehat{m})-F(m)=(\widehat{m}-m)\int_{0}^{1}f\big{(}m+t(\widehat{m}-m)\big{)}dt.$$

Full text

Keywords: nonparametric clustering, kernel estimator, mean shift clustering, plug-in bandwidth, gradient bandwidth.

1 Introduction

Clustering is commonly referred to as the task of finding groups in a set of data points (see [24], [16] or [21]). While intuitively clear, this task is, in fact, far from being accurately defined. The density-based approach attempts to circumscribe this issue by framing the problem into a statistically rigorous setting where the observed data are assumed to be realizations of a random variable, and the clusters are defined with respect to some characteristic of its underlying probability distribution.

In this sense, a clustering procedure should not be limited to simply produce a partition of the observed data; instead, it must allow to obtain a whole-space clustering, that is a partition of the whole sample space [2, 3]. In any case, each methodology is characterized by the way in which the clusters are defined in terms of the true distribution, leading to the concept of ideal population clustering. By serving as a reference “ground truth” to aim at, this concept introduces a benchmark to evaluate the performance of data-based partitions.

The ideal population goal in density-based clustering can be defined in terms of two different paradigms: the model-based approach, where each cluster is associated to a parametric mixture component, and the modal one (see respectively [28] and [30] for some recent reviews). This paper focuses on the latter formulation, whose name stems from the notion of clusters as the “domains of attraction” of the modes of the true density underlying the data [40].

Therefore, in practice density estimation assumes a key role in order to approximate the ideal population goal of modal clustering. While the modal formulation does not preclude using a parametric density estimate as a first step to perform a data-based modal clustering [4, 36], a long-standing practice resorts to nonparametric estimators. Precisely, in this paper the focus lies on those estimators based on kernel smoothing (see e.g. [6] and [45]).

Under- or over-smoothed estimates may lead to deceiving indications about the modal structure of the underlying density function, and this problem is usually quantified through some measure of the discrepancy between the estimate and the target density. In contrast, the aim of this work is to consider nonparametric density estimation as a tool for the final purpose of modal clustering, focusing on an appropriate metric comparing the partitions induced by the true and the estimated distribution.

Our main result provides an asymptotic approximation for the considered metric, which allows to introduce new automatic bandwidth selection procedures specifically designed for nonparametric modal clustering. The accuracy of this approximation and the performance of the new methods in practice, with respect to the proposed error criterion, is extensively studied via simulations, and compared with some plausible competitors.

The rest of the paper is structured as follows. Section 2 formally introduces the modal approach to cluster analysis with reference also to algorithmic details. In Section 3 the distance criterion to target density estimation for modal clustering is presented, along with the main asymptotic result and its consequences. Section 4 contains the setup and results of the numerical experiments. A generalization to the multidimensional setting is discussed in Section 5. Finally, some concluding remarks are stated in Section 6.

2 Background

The connection between groups and density features, established by the modal approach to cluster analysis, allows to characterize the concept of ideal population clustering. Informally, a population cluster can be defined as the domain of attraction of a mode of the density [40]. An attempt to formalize this concept has been done in [3] with the aid of Morse Theory, a branch of differential topology focusing on the large scale structure of an object via the analysis of the critical points of a function (see e.g. [27] for an introduction).

Let us consider a continuous $d$-variate random variable $X$, with probability density function $f\colon\mathbb{R}^{d}\rightarrow\mathbb{R}$. Assume that $f$ is a Morse function, i.e. a smooth enough function having nondegenerate critical points, and denote by $M_{1},\dots,M_{r}$ the modes of $f$ (i.e. its local maxima). For a given initial value $x\in\mathbb{R}^{d}$, an integral curve of the negative density gradient $-\nabla f$ is defined as the path $\nu_{x}:\mathbb{R}\rightarrow\mathbb{R}^{d}$ such that

[TABLE]

The set of points whose integral curve starts at a critical point $x_{0}$ (as $t\to-\infty$) goes under the name of unstable manifold of $x_{0}$ and is defined as

[TABLE]

It has been showed [42] that the class of the unstable manifolds of every critical point of a Morse function yields a partition of the whole space. With these notions at hand, the ideal population clustering $\mathcal{C}=\{\mathcal{C}_{1},\dots,\mathcal{C}_{r}\}$ associated to a density function $f$ is then defined as the set of the unstable manifolds $\{W_{-}^{u}(M_{1}),\dots W_{-}^{u}(M_{r})\}$ of the modes of $f$. By borrowing concepts from terrain analysis, the underlying intuition is that, if $f$ is figured as a mountainous landscape where the modes are the peaks, a modal cluster is the region that would be flooded by a fountain emanating from a peak. When $d=1$, clusters are then unequivocally defined by the locations of the minima points of $f$, which represent the cluster boundaries.

Equivalently, if the integral curves associated to the positive density gradient are considered, then a modal cluster is defined as the set of points whose integral curves converge (as $t\to+\infty$) at the same mode. The concept of modal clusters as the domains of attraction of the density modes stems naturally from this definition. Operationally, a numerical algorithm is needed to find the eventual destination of an initial point, and most of the contributions in this direction take their steps from the mean-shift algorithm [17], essentially a variant of the gradient ascent algorithm. The algorithm transforms an initial point $x^{(0)}$ recursively, and identifies a sequence $(x^{(0)},x^{(1)},x^{(2)},\dots)$ according to an updating mechanism defined as

[TABLE]

where $A$ is a $d\times d$ positive definite matrix chosen to guarantee the convergence to a local maximum of $f$. A partition of the data is therefore obtained by simply grouping together the observations climbing to the same density mode, via mean-shift updates.

From a practical point of view the density $f$ is unknown, therefore an estimate is needed. When working in a nonparametric framework a common choice is given by the kernel density estimator. In the following we focus on the univariate case for ease of exposition and mathematical tractability while the multivariate extension will be addressed in Section 5 below. Let $X_{1},\dots,X_{n}$ be a sample of i.i.d. realizations of $X$. Then, the kernel density estimator is defined by

[TABLE]

where $K$ is the kernel, usually a smooth, non-negative and symmetric function integrating to one, and $h$ is the bandwidth, which controls the smoothness of the density estimate.

While the choice of the function $K$ is known not to have a strong impact in the performance of the estimate [38, Section 3.3.2], choosing $h$ properly turns out to be crucial. A small value of $h$ leads to an undersmoothed density estimate, with the possible appearence of spurious modes, while a too large value results in an oversmoothed density estimate, possibly hiding relevant features.

In order to select the smoothing parameter some measure of the distance between the estimated and the true density is needed. A common choice is the Integrated Squared Error, defined as

[TABLE]

Depending on the observed data, the ISE is itself subject to a random variability that could hinder the problem of bandwidth selection (see [19]). Hence, its expected value

[TABLE]

is alternatively considered as a non-stochastic error distance. The optimal bandwidth $h_{\rm MISE}$ is then defined as $h_{\rm MISE}=\mathop{\rm argmin}_{h>0}{\rm MISE}(h)$.

Since minimization of the MISE does not lead to closed form solutions for the optimal bandwidth, its asymptotic counterpart – the AMISE – is often considered. Both the MISE and the AMISE depend on the true, unknown density function; for this reason several different approaches to estimate them have been proposed. Examples are the ones based on least squares cross validation, biased cross validation or plug-in bandwidth selectors. A comprehensive review of these methods is beyond the scope of this work and, for a complete exposition, readers can refer to [45] or to the more recent book by [6].

3 Density estimation for modal clustering

3.1 Asymptotic bandwidth selection for modal clustering

Bandwidth selectors based on the ISE or akin distances pursue the aim of obtaining an appropriate estimate of the density. However, the goal of modal clustering is markedly different from that of density estimation (see e.g. [12]). In fact, two densities that are close with respect to the ISE may result in quite different clusterings while, on the other hand, densities far away from an ISE point of view could lead to the same partition of the space. A graphical illustration of this idea is provided by Figure 1. The inappropriateness of the ISE, or related distances, depends on its focus on the global characteristics of the density, while modal clustering strongly builds on specific and local features, more closely related to the density gradient or the high-density regions (see also [10]). Therefore, the choice of the amount of smoothing should be tailored specifically for clustering purposes.

So far, the aim of choosing an amount of smoothing for the specific task of highlighting clustering structures has been scarcely pursued in literature. A related idea, although without particular reference to cluster analysis, has been developed by [35], who propose a plug-in type bandwidth selector appropriate for estimation of highest density regions (see also [31] and [14]). Another related work, more focused on the clustering problem, is the one by [15], where the author suggests to consider the self-coverage measure as a criterion for bandwidth selection. Alternatively, the potential adequacy of a bandwidth selected to properly estimate the density gradient has been pointed out informally by [5] and explored numerically by [8]. The theoretical motivation of this suggestion lies on the strong dependence of both the population modal clustering and the mean shift updating mechanism on the density gradient. The suggestion in [9] follows the same rationale and the bandwidth is proposed to be selected as a modification of the normal reference rule for density gradient estimation.

To address the problem of bandwidth selection for modal clustering, an appropriate measure of distance should compare the data-based clustering induced by a kernel density estimate with the ideal population one. Stemming from [3], a natural choice is the distance in measure, where the considered measure here is the probability $\mathbb{P}$ induced by the density $f$. Formally, let $\mathscr{C}=\{C_{1},\dots,C_{r}\}$ and $\mathscr{D}=\{{D}_{1},\dots,{D}_{s}\}$ be two partitions with $r\leq s$ (i.e. possibly different number of groups). The distance in measure between $\mathscr{C}$ and $\mathscr{D}$ is defined as

[TABLE]

where $C\Delta D=(C\cap D^{c})\cup(C^{c}\cap D)$ is the symmetric difference between any two sets $C$ and $D$ and $\mathcal{P}_{s}$ denotes the set of permutations of $\{1,2,\dots,s\}$.

This distance finds an interpretation as the minimal probability mass that would need to be re-labeled to transform one clustering into the other (see Figure 2 for a graphical illustration). In this sense, the second term in (3.1) serves as a penalization for unmatched clusters in one of the clusterings. Practically, this distance conveys the idea that two partitions are similar not when they are physically close, but when the differently-labeled points do not represent a significant portion of the distribution.

It should be noted that the choice of this distance to evaluate the performance of a data-based clustering is not arbitrary. Indeed, many other possibilities are described in [29], but the conclusion of that study is that the distance in measure (called misclassification error there) is “the distance that comes closest to satifying everyone”. Furthermore, in [43] the distance in measure is considered as “the most convenient choice from a theoretical point of view”.

As with the ISE-MISE duality, the distance in measure is a stochastic error distance, so for the purpose of bandwidth selection it seems more convenient to focus on the Expected Distance in Measure

[TABLE]

where $\hat{\mathscr{C}}_{h}$ is the data-based partition induced by $\hat{f}_{h}$ and $\mathscr{C}_{0}$ represents the ideal population clustering. Once the appropriate error distance is defined, the optimal bandwidth $h$ is given by $h_{\rm EDM}=\mathop{\rm argmin}_{h>0}{\rm EDM}(h)$.

As it happened with $h_{\rm MISE}$, it does not seem possible to find an explicit expression for $h_{\rm EDM}$. Hence, our goal will be to obtain an asymptotic form for the EDM that allows to derive a simple approximation to $h_{\rm EDM}$.

To this aim, consider a standard normal random variable $Z$, and denote by $\psi(\mu,\sigma^{2})=\mathbb{E}|\mu+\sigma Z|$ for $\mu\in\mathbb{R}$ and $\sigma>0$. Since $|\mu+\sigma Z|$ has a folded normal distribution [25], it follows that $\psi(\mu,\sigma^{2})$ can be explicitly expressed as

[TABLE]

where $\Phi$ denotes the distribution function of $Z$. This function $\psi$ plays a key role in the asymptotic behavior of the expected distance in measure, as the next result shows (see Appendix A for a proof).

Theorem 1**.**

Assume that $f$ is a bounded Morse function with $r\geq 2$ modes and local minima $m_{1}<\dots<m_{r-1}$, three-times continuously differentiable around each $m_{j}$, that $\int_{-\infty}^{\infty}|x|f(x)dx<\infty$, and that the kernel $K$ is supported on $(-1,1)$, has four bounded derivatives and satisfies $\int_{-\infty}^{\infty}K(x)dx=1$, $\int_{-\infty}^{\infty}xK(x)dx=0$ and $\mu_{2}(K)=\int_{-\infty}^{\infty}x^{2}K(x)dx<\infty$. Define $R(K^{(1)})=\int_{-\infty}^{\infty}K^{(1)}(x)^{2}dx$ and suppose also that $h\equiv h_{n}$ is such that $h\to 0$, $nh^{5}/\log n\to\infty$ and $(nh^{7})^{-1}$ is bounded. Then, ${\rm EDM}(h)$ is asymptotically equivalent to

[TABLE]

where $g^{(k)}$ refers to the $k$-th derivative of a function $g(\cdot)$.

The asymptotically optimal bandwidth $h_{\rm AEDM}$ is then defined as the value of $h>0$ that minimizes ${\rm AEDM}(h)$. Due to the structure of $\psi(\cdot,1)$, minimization of (3.4) is closely related to the problem of minimizing the $L_{1}$ distance in kernel density estimation and, in fact, reasoning as in [20] it is possible to show that $h_{\rm AEDM}$ is of order $n^{-1/7}$. Unfortunately, as it happened with $h_{\rm EDM}$, it seems that neither $h_{\rm AEDM}$ admits an explicit representation hence, to get further insight into the problem of optimal bandwidth selection for density clustering, it appears necessary to rely on a tight upper bound for ${\rm AEDM}(h)$.

To find such a bound it is useful to note that many properties of $\psi(u,1)$ are given in [13, Ch. 5], and can be translated to our function of interest by taking into account that $\psi(\mu,\sigma^{2})=\sigma\psi(\mu/\sigma,1)$. It follows that $\psi(\mu,\sigma^{2})$ is symmetric with respect to $\mu$, nondecreasing for $\mu>0$ and convex, attaining its minimum at $\mu=0$ so that $\psi(\mu,\sigma^{2})\geq\psi(0,\sigma^{2})=(2/\pi)^{1/2}\sigma$ for all $\mu\in\mathbb{R},\sigma>0$.

By taking into account that $e^{-\mu^{2}/(2\sigma^{2})}$ and $1-2\Phi(-\mu/\sigma)$ are both bounded by 1, [13] also noted that

[TABLE]

for all $\mu\in\mathbb{R},\sigma>0$. However, a tighter bound for small values of $\mu$ is given in the next lemma.

Lemma 1**.**

The bound $\psi(\mu,\sigma^{2})\leq(2/\pi)^{1/2}\sigma+(2\pi)^{-1/2}\mu^{2}/\sigma$ holds for all $\mu\in\mathbb{R}$ and $\sigma>0$.

The bound in Lemma 1 is tighter than (3.5) whenever $|\mu|\leq(2\pi)^{1/2}\sigma$, but the situation reverses for bigger values of $|\mu|$, so that none of the two bounds is uniformly better (see Figure 3) hence we should keep track of both of them. They lead to upper bounds for the asymptotic EDM.

Corollary 1**.**

Under the conditions of Theorem 1, the asymptotic EDM satisfies ${\rm AEDM}(h)\leq\min\{{\rm AB}1(h),{\rm AB}2(h)\}$ for all $h>0$, where

[TABLE]

Here, $b=\sum_{j=1}^{r-1}b_{j}$ and $a_{\ell}=\sum_{j=1}^{r-1}a_{j\ell}$ and for $\ell=1,2$, where

[TABLE]

The minimizers of ${\rm AB}1(h)$ and ${\rm AB}2(h)$ can be computed explicitly, and are given by

[TABLE]

3.2 Some remarks

In this section we discuss in more depth some of the results derived in Section 3.1. The aim is to provide insights on the behavior of the approximations and bandwidth selectors and to discuss possible competitors.

Remark 1*.*

Theorem 1 provides an asymptotic expression for the EDM that is valid as long as the true density has two or more modes. When the true density is unimodal ($r=1$), expression (3.4) is not well-defined. However, under the assumptions of the theorem the kernel estimator is also unimodal with probability one for big enough $n$. Thus, asymptotically the distance in measure would be identically zero, hence the AEDM formula would remain valid under the usual convention setting $\sum_{j=1}^{0}=0$.

Moreover, for unimodal densities the numerical work in Section 4 suggests that there exists $h_{0}>0$ such that ${\rm EDM}(h)=0$ for all $h\geq h_{0}$. Hence, in that case it seems sensible to define $h_{\rm EDM}=\inf\{h>0\colon{\rm EDM}(h)=0\}$.

Remark 2*.*

A natural estimator of the density first derivative is the first derivative of the kernel density estimator. For this estimator it is possible to define the MISE as in (2.1), and to consider its minimizer $h_{\rm MISE,1}$ and its asymptotic approximation $h_{\rm AMISE,1}$ (see [39] and [7]). The bandwidths (3.6) and (3.7) share the same order as $h_{\rm AMISE,1}$, whose expression is given by

[TABLE]

with $R(f^{(3)})=\int_{-\infty}^{\infty}f^{(3)}(x)^{2}dx$. This consideration strengthens the intuition, outlined in Section 3.1, that (3.8) could be an adequate bandwidth choice for modal clustering purposes.

Remark 3*.*

By explicitly plugging expression (3.1) for $\psi$ into (3.4), it is easily seen that the AEDM can be decomposed into two summands. Studying their behavior, as a function of $h$, it can be checked that the first term decreases when $h\to 0$ while the second one tends to increase, and viceversa for $h$ taking large values. A similar trade-off occurs with the decomposition of the AMISE into the Asymptotic Integrated Squared Bias and the Asymptotic Integrated Variance, which are minimized for diverging values of $h$.

Remark 4*.*

If the true density is exactly symmetric around its minimum, the considerations in the previous item do not hold anymore. Symmetry around a minimum $m$ implies $f^{(k)}(m)=0$, for any odd value of $k$. Therefore the first summand of the AEDM expression, related to the bias, vanishes, leading to a monotonically decreasing behavior of the AEDM itself. A similar anomaly was observed in the related problem of mode estimation in [11]: if the true density is symmetric around its mode, then Chernoff’s mode estimator is unbiased. Hence, in some special cases symmetry plays a certain role in the performance of these smoothing methodologies.

Remark 5*.*

The derived bandwidths depend on some unknown quantities such as the true density, its local minima and its second and third derivatives. In order to be of practical use we shall resort to plug-in strategies, that is, data-based bandwidth selectors will be proposed in the next section by substituting the aforementioned unknown quantities with pilot estimates. This is the same procedure that is commonly adopted when considering the plug-in bandwidth selector $\hat{h}_{\rm PI,1}$ for density gradient estimation (see [23] and [5]).

It should be noted that this allows, from a practical point of view, to circumvent the issue about the perfect symmetry around a minimum since, by using a nonparametric pilot estimate of the third derivative, it is highly unlikely to encounter a similar situation in practice.

4 Numerical results

The idea of estimating the density for clustering purposes, via the minimization of the expected distance in measure – or its asymptotic counterpart – is explored in this section via simulations. All the analyses have been performed in the R environment [32] with the aid of the ks [41], meanShiftR [26], clue [22], and multimode [1] packages.

A total of $B=1000$ samples for each of the sizes $n\in\{100,1000,10000\}$ are generated from the univariate densities depicted in Figure 4 and whose parameters are reported in Appendix B. The selected densities are designed to illustrate different modal structures to encompass different possible behaviors from a clustering perspective.

The first goal of the study was to evaluate the quality of the asymptotic approximation of the EDM and the behavior of the two bounds derived in Corollary 1. Since an explicit expression for the EDM was not available, we obtained a Monte Carlo approximation based on the $B=1000$ synthetic samples.

The plots displayed in Tables 1 to 5 show the behavior of the asymptotic approximations, with respect to the EDM, as a function of the bandwidth $h$. As expected, the approximations improve as the sample size increases. The two bounds show a quite different behavior, with characteristics that reflect the theoretical properties pointed out in Section 3.2. The first bound is closer to the AEDM in uniform terms, but despite having a diverging behavior for large $h$ the second bound is usually closer to the AEDM around the location of the minimizer $h_{\rm AEDM}$.

With regard to the EDM, it presents a nearly flat pattern around its minimizer, thus suggesting a range of plausible bandwidths with very similar performance as the optimal one. This is especially true for densities with a simpler modal structure, captured by the kernel estimate for a wide range of bandwidth values, which suggests that bandwidth selection for modal clustering is, in fact, easier than for density estimation.

To appreciate how much is lost by changing the target from the optimal $h_{\rm EDM}$ to the oracle surrogates $h_{\rm AEDM}$ and $h_{\rm MISE,1}$, the first three lines in each table also present the values for the corresponding EDM, all computed under a full knowledge of the density and its involved features. By construction, ${\rm EDM}(h_{\rm EDM})$ is the lowest of these values and, being derived as an asymptotic approximation, the oracle $h_{\rm AEDM}$ stands close to this optimal value, especially for larger sample sizes. However, it is remarkable that $h_{\rm MISE,1}$, despite being based on a different optimality criterion, also leads to comparable or even improved results over $h_{\rm AEDM}$ in terms of the EDM.

As a second goal, we propose new data-based bandwidth selectors specifically designed for modal clustering purposes. The first step consists in estimating the number of local minima, and their location. This is achieved by numerically finding the roots of a pilot estimate of $f^{(1)}$, constructed as the derivative of the kernel density estimator using the plug-in gradient bandwidth $\hat{h}_{\rm PI,1}$. Then, similarly, using plug-in bandwidths for $f$, $f^{(2)}$ and $f^{(3)}$, we obtain pilot estimates of these involved unknown functions at the estimated local minima. These quantities are subsequently plugged-in in the formulas of the AEDM, AB1 and AB2, and the minimizers of the resulting estimated criteria are found; in the case of the estimated AEDM by numerical minimization, and according to expressions (3.6) and (3.7) for AB1 and AB2 respectively. The data-based bandwidths thus obtained are denoted $\hat{h}_{\rm AEDM}$, $\hat{h}_{\rm AB1}$ and $\hat{h}_{\rm AB2}$, respectively.

Occasionally (although rarely) the first step in the procedure above yielded a single mode, and then the AEDM was undefined. In those cases, and according to the rationale exposed in Remark 1, a sensible choice for $h$ is the critical bandwdith proposed by [37],

[TABLE]

so in that case we set $\hat{h}_{\rm AEDM}=\hat{h}_{\rm AB1}=\hat{h}_{\rm AB2}=\hat{h}_{\rm crit}$.

Tables 1 to 5 also contain the Monte Carlo averages and standard deviations of the distances in measure obtained when performing modal clustering using the bandwidth selectors $\hat{h}_{\rm AEDM}$, $\hat{h}_{\rm AB1}$ and $\hat{h}_{\rm AB2}$. For completeness, their performance is also compared to that of $\hat{h}_{\rm PI,1}$, which so far probably represents their most sensible competitor in the clustering framework (see [8]).

In general, $\hat{h}_{\rm AB1}$ and $\hat{h}_{\rm AB2}$ led to more accurate clusterings than $\hat{h}_{\rm AEDM}$, with a slight preference for $\hat{h}_{\rm AB1}$. The gradient-based bandwidth $\hat{h}_{\rm PI,1}$, in turn, not only produces competitive results, but its Monte Carlo average distance in measure appears lower than the one produced by the asymptotic EDM minimizers. In fact, a deeper insight into the standard errors of the obtained distances shows that $\hat{h}_{\rm AEDM}$, as well as $\hat{h}_{\rm AB1}$ and $\hat{h}_{\rm AB2}$, produce more variable results. The higher variability seems to be due to the sensitivity of the minimizers to the plugged in pilot estimates, which strongly depend on local features of the density. Some further investigations, not fully reported here, suggest that the main responsible for this behaviour is not the pilot estimate of the local minima but the pilot density derivatives estimates at the minimum points. Also, due to the use of different pilot bandwidths to estimate the unknown $m_{j}$, $f^{(2)}$, and $f^{(3)}$, it may occur, indeed, that $\hat{f}^{(2)}(\hat{m}_{j})$ assumes even negative values. On the other hand, while relying as well on some plug-in estimates, the gradient-based bandwidth $\hat{h}_{\rm PI,1}$ produces more robust clusterings, as the quantities to be estimated refer conversely to global features of the density. As expected, this diverging behavior tends to vanish with increasing sample size since the asymptotic approximations improve. As a confirmation, with $n=10000$, all the considered bandwidths perform comparably.

5 Multidimensional generalization

The concepts discussed so far refer to the one-dimensional setting where a mathematically rigorous treatment is feasible. The multidimensional generalization poses some difficulties since obtaining an asymptotic approximation of the EDM appears far from trivial. Hence, in order to gain some insight into the problem of selecting the amount of smoothing for nonparametric clustering in more than one dimension, some numerical comparisons are performed assuming the true density as known.

Denote by $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ the true density function and by

[TABLE]

its kernel estimate based on a sample $\mathbf{X}_{1},\dots,\mathbf{X}_{n}$ and indexed by a symmetric positive definite $d\times d$ bandwidth matrix $\mathbf{H}$. The problem of bandwidth selection is considered by studying the EDM between the clustering induced by the kernel estimate $\hat{\mathscr{C}}_{\mathbf{H}}$ and the ideal population clustering $\mathscr{C}_{0}$. These clusterings are not so easily identifiable as in the unidimensional setting, due to the arbitrary forms that the cluster boundaries may adopt, however an approximation of the distance in measure $d(\hat{\mathscr{C}}_{\mathbf{H}},\mathscr{C}_{0})$ can be computed by resorting to a discretization scheme as follows (see [8] for further details):

Take a grid over the sample space and rule the grid by considering hyper-rectangles centered at each grid point. 2. 2.

Assign a cluster membership to each grid point by running a population version of the mean-shift algorithm i.e. using the true density. This produces a discretized version of $\mathscr{C}_{0}$. 3. 3.

Similarly, obtain the data-based partition $\hat{\mathscr{C}}_{\mathbf{H}}$ induced by $\hat{f}_{\mathbf{H}}.$ 4. 4.

Compute the probability mass of each single hyper-rectangle in $\mathscr{C}_{0}$. 5. 5.

Compute the distance in measure as in (3.2) where the involved probabilities are evaluated based on the previous step.

For the multidimensional simulation study, a total of $B=1000$ samples for each of the sizes $n\in\{100,1000\}$ were generated from the bivariate densities whose contour plots are shown in Figure 5 and described in Appendix B. The densities have been chosen to generalize the settings M1 and M5 included in the univariate study.

Three different parametrizations for the bandwidth matrix were considered: a scalar bandwidth $\mathbf{H}=h^{2}\mathbf{I},$ with $\mathbf{I}$ the identity matrix, a diagonal bandwidth $\mathbf{H}=\text{diag}(h_{1}^{2},h_{2}^{2})$, and a full, unconstrained bandwidth matrix $\mathbf{H}.$ For density and density derivative estimation, [44] and [7] showed that the use of the simplest scalar bandwidth can be quite detrimental in practice, a diagonal bandwidth may suffice in some scenarios, but in general it is advantageous to employ unconstrained bandwidth matrices (see also [6]). However, such results have never been obtained in a modal clustering framework; thus one of the goals of this simulation study is to examine how the bandwidth matrix parametrization affects the performances of the procedures.

Using the synthetic samples from each density in the study, it was possible to obtain a Monte Carlo estimate of the (discretized version of the) EDM, which was then minimized over the class of scalar, diagonal and unconstrained bandwidth matrices. The EDM was computed also for the MISE-optimal bandwidth for density gradient estimation over the same matrix classes. In both cases, the true density as well as all the involved quantities were assumed to be known. The EDM minimizers were determined numerically, by running the procedure over a grid of sensible values of the entries, while the optimal matrices for gradient estimation were determined as in [7].

The results are reported in Tables 6 and 7. Clustering based on the optimal bandwidth according to the EDM is very accurate in both of the considered examples, and improves considerably for increasing sample size. The use of more complex bandwidth parametrizations does not seem worth for modal clustering since results obtained with a full, unconstrained bandwidth matrix are comparable with those obtained with a scalar bandwidth, while the latter requires a substantially smaller computational effort.

In the multidimensional setting, the gradient bandwidth is quite competitive in terms of EDM, as in the univariate case. Again the comparable performance of unconstrained bandwidth matrices does not seem to justify the use of more complex parametrizations.

6 Conclusions

The modal clustering methodology provides a framework to perform cluster analysis with a clear and explicit population goal. It allows clusters of arbitrary shape and size, which can be captured by means of a nonparametric density estimator. In this context, the distance in measure represents a natural and easily interpretable error criterion. Therefore, in this paper we have presented an asymptotic study of this criterion for the case where density estimates of kernel type are employed to obtain a whole-space clustering via the mean shift algorithm.

Our asymptotic approximations are useful to gain insight into the fundamental problem of bandwidth selection for modal clustering and, at the same time, serve as the basis to propose practical data-based bandwidth choices specifically designed for clustering purposes.

The finite-sample performance of the new proposals was investigated in a thorough simulation study, and compared to the oracle bandwidths i.e. the optimal choices when the true population is fully known. The gradient bandwidth, designed for the closely related problem of density gradient estimation, was also included as a natural competitor in the study.

The results of this simulation study have suggested that all the methods perform quite satisfactorily, and exhibit a very similar behavior for large sample sizes. For smaller samples, the performance of the gradient bandwidth was rather remarkable, since it obtained comparable or even better results than the new proposals, even without being specifically conceived for modal clustering.

This phenomenon resembles the conclusions obtained in [34] regarding the related problem of level set estimation. There, it was shown that the traditional bandwidth selectors for density estimation often outperformed more sophisticated methods designed for level set estimation purposes. The common pattern in both situations is that the optimal choices for the specific problems (level set estimation and modal clustering, respectively) depend on very subtle local features of the unknown density function, which are difficult to estimate, so that choices based on a more global, yet somehow related, perspective represent a sensible alternative.

Appendix A Proofs

Proof of Theorem 1.

From Theorem 4.1 in [3] it follows that, with probability one, there exists $n_{0}\in\mathbb{N}$ such that the kernel density estimator $\widehat{f}_{h}$ has the same number of local minima as $f$ for all $n\geq n_{0}$. Let us denote by $\widehat{m}_{h,1}<\dots<\widehat{m}_{h,r-1}$ the local minima of $\widehat{f}_{h}$. Then, the expected distance in measure between the data-based clustering $\widehat{\mathscr{C}}_{h}$ and the population clustering $\mathscr{C}_{0}$ can be written as

[TABLE]

Write, generically, $\widehat{m}$ and $m$ for any of the estimated and true local minima. A Taylor expansion with integral remainder allows to write

[TABLE]

The assumptions imply that $\widehat{m}\to m$ almost surely [see, for instance, 33] and, since $f$ is bounded and continuous, this readily yields \int_{0}^{1}f\big{(}m+t(\widehat{m}-m)\big{)}dt\to f(m) almost surely, which entails that $\mathbb{E}|F(\widehat{m})-F(m)|\sim f(m)\mathbb{E}|\widehat{m}-m|$. The result then follows from Equation (2.6) in [18], where the asymptotic form of $\mathbb{E}|\widehat{m}-m|$ is given. ∎

Proof of Lemma 1.

From $\psi(\mu,\sigma^{2})=\sigma\psi(\mu/\sigma,1)$, it suffices to show that $\psi(u,1)\leq(2/\pi)^{1/2}+(2\pi)^{-1/2}u^{2}$ for $u\geq 0$. From the definition of $\psi$, this is equivalent to proving that $\alpha(u)\leq 1$, where $\alpha(u)=e^{-u^{2}/2}+u\int_{0}^{u}e^{-z^{2}/2}dz-u^{2}/2$. Since $\alpha(0)=1$, it is enough to show that $\alpha$ is nonincreasing, but this immediately follows from the fact that $\alpha^{\prime}(u)=\int_{0}^{u}e^{-z^{2}/2}dz-u$. ∎

Appendix B Parameter settings

In the following the parameter settings of the densities selected for the simulations are presented. Since all the densities are mixture of Gaussian models, we adopt the usual notation where, for a given $k$ component, $\pi_{k}$ represent the k-th mixture weight, $\mu_{k}$ and $\sigma_{k}^{2}$ ($\Sigma_{k}$ for the bivariate models) the mean and variance (covariance matrix).

B.1 Unidimensional parameter settings

B.1.1 Density M1

B.1.2 Density M2

B.1.3 Density M3

B.1.4 Density M4

B.1.5 Density M5

B.2 Bidimensional settings

B.2.1 Asymmetric bimodal

B.2.2 Trimodal