Density Level Sets: Asymptotics, Inference, and Visualization

Introduction

Estimating the level sets of a probability density function has a wide range of applications, including anomaly detection (outlier detection) (Breunig et al., 2000; Hodge and Austin, 2004), two-sample comparison (Duong et al., 2009), binary classification (Mammen and Tsybakov, 1999), and clustering (Rinaldo and Wasserman, 2010; Rinaldo et al., 2012). In this paper, we study the problem of estimating the level set

where $p_{h}$ is the expected kernel density estimator with bandwidth $h$, a smoothed version of the underlying density $p$. Using $D_{h}$ (and thus $p_{h}$) has several advantages, which we discuss in detail in Section 2.2. Figures 5 and 10 illustrate the kind of confidence sets and visualizations we will develop in this paper.

A commonly used estimator of the density level-set is the plug-in estimator $\widehat{D}_{h}=\left\{{x:\;\widehat{p}_{h}(x)=\lambda}\right\}$, where $\widehat{p}_{h}$ is the kernel density estimator or some other density estimator. There is a large literature for level sets (and upper level sets, which replace $=\lambda$ with $\geq\lambda$) that focuses on the consistency, rates of convergence (Polonik, 1995; Tsybakov, 1997; Walther, 1997; Cadre, 2006; Cuevas et al., 2006) and minimaxity (Singh et al., 2009) of such estimators under various error loss functions.

Recent results on statistical inference for level sets include Jankowski and Stanberry (2012) and Mammen and Polonik (2013). Statistical inference is challenging in this setting because the estimand is a set and the estimator is a random set (Molchanov, 2005). Mason and Polonik (2009) establish asymptotic normality for upper level sets when the loss function is the measure of the set difference. However, it is unclear how to derive a confidence set from this result.

Another challenge of level set estimation is that we cannot directly visualize the level sets when the dimension of the data $d$ is larger than 3. One approach is to construct a level-set tree, which shows how the connected components for the upper level sets bifurcate when we gradually increase $\lambda$ (Stuetzle, 2003; Klemelä, 2004, 2006, 2009; Stuetzle and Nugent, 2010; Kent et al., 2013). The level-set tree reveals topological information about the level sets but loses geometric information.

In this paper, we propose solutions to all of these problems. Our main contributions can be summarized as follows.

We derive the limiting distribution of ${\sf Haus}(\widehat{D}_{h},D_{h})$ (Theorem 3).

We develop two bootstrap-based methods to construct confidence regions for $D_{h}$ (Section 4).

We prove that both bootstrap methods are valid (Theorem 28 and 5).

We devise a visualization technique that preserves the geometric information for density level sets (Section 5).

Related Work. Early work on density level set focuses on proving the consistency or the rate of convergence under various metrics. See e.g. Polonik (1995); Tsybakov (1997); Walther (1997); Cuevas et al. (2006); Rinaldo and Wasserman (2010). However, none of these derives a limiting distribution for the density level sets. To our knowledge, the only paper that considers limiting distributions is Mason and Polonik (2009), proving asymptotic normality under a generalized integrated distance. However, this metric cannot be used to construct a confidence set for density level sets since the asymptotic distribution involves the true density, which is unknown. Estimating the level set is also related to support estimation, see e.g. Cuevas and Rodríguez-Casal (2004) and Cuevas (2009).

Jankowski and Stanberry (2012) and Mammen and Polonik (2013), both provide methods for constructing confidence sets for the density level sets using the variation of the density function. Our approach is similar to theirs but is based on Hausdorff distance. We will compare their methods to ours in Section 4.2.

Outline. We begin with a short introduction to density level sets along with some useful geometric concepts in Section 2. In Section 3, we derive the limiting distribution of the Hausdorff distance between the estimated and true level sets. In Section 4, we construct a valid confidence set for density level sets. In Section 5, we devise a visualization method for density level sets that is simple to interpret and efficient to compute. (We provide an R package that implements our visualization method.) We summarize our results and discuss related problems in Section 6.

Technical Background

Let $X_{1},\cdots,X_{n}$ be a random sample from an unknown, continuous density $p(x)$. We define the density level set by

for some $\lambda>0$. Note that in the literature the term level sets is sometimes used for the set $\left\{{x:\;p(x)\geq\lambda}\right\}$; we call the latter the upper level set to distinguish it from the level set in equation 2. Thus, the level sets (in our terminology) are the boundaries to the upper level sets, under mild smoothness assumptions (e.g., assumption G below).

We assume that $\lambda$ is a fixed, positive value. A plug-in estimate for $D$ is

where $\widehat{p}_{h}$ is the kernel density estimator (KDE),

2 Smoothed Density Level Set

In this paper, we focus on inference for the level sets of a smoothed version of $p$, specifically:

where $K_{h}(x)=\frac{1}{h^{d}}K\left(\frac{\|x\|}{h}\right)$ and $\star$ denotes convolution. We denote the $\lambda$ level set of $p_{h}$ by

Note that although we focus on estimating $D_{h}$, we allow $h=h_{n}\rightarrow 0$ as $n\rightarrow\infty$.

Here, we argue that $p_{h}$ (and hence $D_{h}$) is a better target for level-set estimation than $p$ (and hence $D$). For simplicity, in this section we focus on the upper level sets $L_{h}(\lambda)=\{x:\ p_{h}(x)\geq\lambda\}$ and $L(\lambda)=\{x:p_{h}(x)\geq\lambda\}$, but the gist of the argument remains the same in either case.

Our arguments are: (i) $p_{h}$ always exists while $p$ may not even exist; (ii) $\widehat{p}_{h}$ can be naturally viewed as an estimator for $p_{h}$; (iii) $p_{h}$ has all the salient structure of $p$ but is easier to estimate than $p$; (iv) While the bias $p_{h}(x)-p(x)$ may be analyzed theoretically, in practice, it cannot be accurately estimated. It is better to just be clear that we are really estimating $p_{h}$.

As long as $h\geq(\log n/n)^{1/d}$, $P^{n}(\sup_{x}|\widehat{p}_{h}(x)-p_{h}(x)|\geq\epsilon)\to 0$, and hence we can uniformly consistently estimate $p_{h}$, whether we keep $h$ fixed or let it tend to 0. The same is not true for $p$. In fact, $P$ may not even have a density $p$. The bias cannot be uniformly estimated in a distribution-free way.

Regarding (iii): The left plot in Figure 1 shows a density $p$. The blue points at the bottom show the upper level set $L=\{x:\ p\geq 0.05\}$. The right plot shows $p_{h}$ for $h=0.2$ and the blue points at the bottom show the upper level set $L_{h}=\{x:\ p_{h}\geq 0.05\}$. The smoothed out density $p_{h}$ is biased and the upper level set $L_{h}$ loses the small details of $L$. But these small details are the least estimable, and $L_{h}$ captures the principal structure of $L$. In addition, when $L$ is smooth (having positive $\mu$-reach; Chazal and Lieutier 2005; Chazal et al. 2009) and $h$ is sufficiently small, $L_{h}$ and $L$ will be topologically similar, in the sense that a small expansion of $L_{h}$ is homotopic to $L$ (Chazal and Lieutier 2005; Chazal et al. 2009; Genovese et al. 2014). The $\mu$-reach and the concept of being nearly homotopic can be found in section 3.2 of Genovese et al. (2014) and section 4 of Chazal et al. (2009).

As a second example, let $P=(1/3)\phi(x;-5,1)+(1/3)\delta_{0}+(1/3)\phi(x;5,1)$ where $\phi$ is a Normal density and $\delta_{0}$ is a point mass at 0. Of course, this distribution does not even have a density. The left plot in Figure (2) shows the density of the absolutely continuous part of $P$ with a vertical line to who the point mass. The right plot shows $p_{h}$, which is a smooth, well-defined density. Again the blue points show the level sets. As before $p_{h}$ and $L_{h}$ are slightly biased, but they are also well-defined. And $\widehat{p}_{h}$ and $\widehat{L}_{h}$ are accurate estimators of $p_{h}$ and $L_{h}$, respectively. Moreover, $L_{h}$ captures the most important qualitative information about $L$, namely, that there are three connected components, one of which is small. These examples show that $p_{h}$ — and hence $D_{h}$ — is a sensible target of inference.

Lastly, we would like to point out that the idea of viewing $p_{h}$ as the estimand is not new. The “scale space” approach to smoothing explicitly argues that we should view $\widehat{p}_{h}$ as an estimate of $p_{h}$, and $p_{h}$ is then regarded as a view of $p$ at a particular resolution. This idea is discussed in detail in Chaudhuri et al. (2000); Chaudhuri and Marron (1999); Godtliebsen et al. (2002).

If one really wants to focus on making inference for $D$, the level set of the original density, then we need to undersmoothUndersmoothing the density estimate to make statistical inferences is a common practice in nonparametric statistics; see e.g. page 89 of Wasserman (2006). so that the bias will not affect the limiting distribution. This leads to estimates of $D$ that are highly variable. We believe that an accurate confidence set for $D_{h}$ is more useful than a poor confidence set for $D$.

These arguments explain why we regard $p_{h}$ rather than $p$ as the estimand. But these arguments do not tell us how to choose $h$. Bandwidth selection is always a challenge and in this paper we mainly use Silverman’s rule of thumb (Silverman (1986)).

3 Geometric Concepts

Let $\pi_{A}(x)$ be the projection of a point $x$ onto a set $A$. Note that $\pi_{A}(x)$ may not be unique. The distance induced by the projection is

where $A\oplus\epsilon=\bigcup_{x\in A}B(x,\epsilon)$ and $B(x,\epsilon)=\{y:\;\|x-y\|\leq\epsilon\}$. The Hausdorff distance is a generalized version of the ${\cal L}_{\infty}$ metric for sets.

The reach of a set $M$ (Federer 1959; Cuevas 2009, also known as condition number Niyogi et al. 2008 or minimal feature size Chazal and Lieutier 2005) is the largest distance from $M$ such that every point within this distance to $M$ has a unique projection onto $M$. i.e.

Note that $\pi_{A}(x)$ is unique if $0<d(x,A)\leq\operatorname{{\sf reach}}(A)$. Another way to understand the reach is as the largest radius of a ball that can freely move along $M$; see Figure 3 for an example. In some cases, the reach is the same as the smallest radius of curvature on $M$. The reach plays a key role in relating the Hausdorff distance to the empirical process. Note that the reach is closely related to ‘rolling properties’ and ‘$\alpha$-convexity’; see Cuevas (2009), Cuevas et al. (2012) and appendix A of Pateiro-López (2008).

Finally, two smooth sets $A$ and $B$ are called normal compatible (Chazal et al., 2007) if the projection between $A$ and $B$ are one to one and onto; see Figure 4 for an example. When $A$ and $B$ are normal compatible, the Hausdorff distance has the simpler form

Asymptotic Theory

In this section, we derive the asymptotic theory for ${\sf Haus}(\widehat{D}_{h},D_{h})$. Note first that for two sets $A$ and $B$, the Hausdorff distance satisfies the inclusion property

From this, it follows that when we have a set estimator $\widehat{A}_{n}$ and a parameter of interest $A$, then for any $\alpha>0$, the set

is a $1-\alpha$ confidence set for $A$, where $\operatorname{{\sf Quantile}}(X,\alpha)$ is the $\alpha$-quantile of random variable $X$. Thus, whenever we can approximate the distribution of ${\sf Haus}(\widehat{A}_{n},A)$, we can construct a confidence set for $A$. Note that Chen et al. (2014b) and Chen et al. (2014a) have used this property to construct confidence sets, but neither paper used this property to full effect.

We define the sup norm using derivatives up to $r$-th order by:

Let $\alpha=(\alpha_{1},\cdots,\alpha_{d})$ be an multi-index such that each $\alpha_{j}$ is a non-negative integer and $|\alpha|=\alpha_{1}+\cdots+\alpha_{d}$. We define

We now state are main assumptions, for an arbitrary density $q$. When we use the assumptions in what follows, we will take $q$ to be $p_{h}$.

The kernel function $K\in\mathbf{BC}^{3}$ and is symmetric, non-negative, and

for all multi-indices $\alpha$ satisfying $|\alpha|\leq 3$.

The kernel function $K$ and its partial derivative satisfies condition $K_{1}$ in Giné and Guillou (2002). Specifically, let

Assumption (G) appears in Molchanov (1991); Tsybakov (1997); Walther (1997); Molchanov (1998); Cadre (2006); Mammen and Polonik (2013); Laloe and Servien (2013). For a smooth density $q$, (G) holds whenever the specified level $\lambda$ does not coincide with the density value for a critical point.

Assumption (K1) is to guarantee that the variance of the KDE is bounded and to ensure that $p_{h}\in\mathbf{BC}^{3}$. This assumption is very common in statistical literature, see e.g. Wasserman (2006). Assumption (K2) is to regularize the complexity of the kernel function so that the supremum norm for kernel functions and their derivatives can be bounded in probability. Similar assumption appears in Einmahl and Mason (2005) and Genovese et al. (2014). The Gaussian kernel and many compactly supported kernels satisfy both assumptions.

An immediate result from assumption (G) is the smoothness of the density level set. This smoothness property will be used to understand the distribution of ${\sf Haus}(\widehat{D}_{h},D_{h})$.

Assume a density $p_{h}\in\mathbf{BC}^{2}$ satisfies condition (G) and let $D_{h}$ denote the level set for $p_{h}$ at $\lambda$. Then

Moreover, let $q\in\mathbf{BC}^{3}$ be another density function and define $D(q)$ as the level set for $q$ at level $\lambda$. When $\|p_{h}-q\|^{*}_{2,\max}$ is sufficiently small,

$\operatorname{{\sf reach}}(D(q))=\min\left\{\frac{\delta_{0}}{2},\frac{g_{0}}{\|p_{h}\|^{*}_{2,\max}}\right\}+O(\|p_{h}-q\|^{*}_{2,\max})$.

$D_{h}$ and $D(q)$ are normal compatible.

The proof is given in the supplementary materials. Lemma 1 is very similar to Theorem 1 and 2 in Walther (1997). Essentially, this lemma shows that the level set $D$ is smooth and that whenever two smooth densities are sufficiently close, their level sets will both be smooth, be close to each other, and have one-to-one and onto normal projections between them.

Assume (K1–K2) and (G) hold for $p_{h}$. Let $D_{h}$ and $\widehat{D}_{h}$ be the density level sets with level $\lambda$ for $p_{h}$ and $\widehat{p}_{h}$. Define the function

with $x\in D_{h}$. If $\frac{\log n}{nh^{d+2}}\rightarrow 0,\,h\rightarrow 0$, then

A key element for the proof of Lemma 19 is the smoothness of $D_{h}$ and $\widehat{D}_{h}$, which relies on Lemma 1. This smoothness allows us to approximate the local difference by an empirical process.

Lemma 19 shows that the projected distance to the level sets can be approximated by a stochastic process (the empirical process) defined on a smooth manifold. Specifically, Lemma 19 shows that the projection distance can be approximated by an empirical process on certain functions $f_{x}$, where $x\in D_{h}$. The level sets $D_{h}$ now acts as an index set. Thus, we define the function space

Assume (K1–K2) and (G) holds for $p_{h}$. Let $D_{h}$ and $\widehat{D}_{h}$ be the density level sets with level $\lambda$ for $p_{h}$ and $\widehat{p}_{h}$. Then when $\frac{\log n}{nh^{d+2}}\rightarrow 0,\,h\rightarrow 0$, the Hausdorff distance satisfies

The proof of Theorem 3 depends on two geometric observations. First, the empirical approximation in Lemma 19, shows that the local difference is approximately the same as an empirical process, and hence the maximum local difference is approximated by the maximum of the empirical process. Second, the normal compatibility between $\widehat{D}_{h}$ and $D_{h}$ guaranteed by Lemma 1, which implies that maximal of local difference is the same as the Hausdorff distance.

Theorem 3 shows that the Hausdorff distance ${\sf Haus}(\widehat{D}_{h},D_{h})$ can be approximated by a maximum over a certain Gaussian process. Note that we cannot directly use this theorem to construct a confidence set for $D_{h}$ since the Gaussian process is defined on $D_{h}$, which is unknown. Later we will use the bootstrap to approximate this limiting distribution and construct a confidence set.

Statistical Inference

We now show that we can construct valid confidence sets for $D_{h}$ with the bootstrap. A set $S_{n,1-\alpha}$ is called an asymptotically valid confidence set for $D_{h}$ if

where $r_{n}\rightarrow 0$ as $n\rightarrow\infty$. We propose two methods for constructing a confidence set, and we will show that they are both asymptotically valid.

The first approach is to use the Hausdorff distance between the level sets. Let $W_{n}={\sf Haus}(\widehat{D}_{h},D_{h})$ and define

where $F_{A}$ denotes the cdf for a random variable $A$. Then, it is easy to see that

We use the bootstrap to estimate $w_{1-\alpha}$.

Let $X^{*}_{1},\cdots,X^{*}_{n}$ be a bootstrap sample from the empirical measure. Let $\widehat{p}_{h}^{*}$ denote the KDE using the bootstrap sample, and $\widehat{D}^{*}_{n}$ the corresponding level set. We define $W^{*}_{n}={\sf Haus}(\widehat{D}^{*}_{n},\widehat{D}_{h})$ and

Then the bootstrap confidence set is $\widehat{D}_{h}\oplus\widehat{w}_{1-\alpha}$.

Assume (K1–K2) and (G) holds for $p_{h}$ and $\frac{\log n}{nh^{d+2}}\rightarrow 0,\,h\rightarrow 0$. Let $D_{h}$ and $\widehat{D}_{h}$ and $\widehat{D}^{*}_{n}$ be the density level set with level $\lambda$ for $p_{h}$ and $\widehat{p}_{h}$ and $\widehat{p}_{h}^{*}$. Then there exist $\mathcal{X}_{n}$ such that

An intuitive explanation for Theorem 28 is that as $n$ goes to infinity, the bootstrap process converges to the same Gaussian process as the empirical process – thus, they share the same Berry-Esseen bound.

2 Method 2: Supremum Loss

The second approach is to use the supremum norm of the KDE and impose an upper and lower bound around the density level.

Assume (K1–K2) and (G) holds for $p_{h}$ and $\frac{\log n}{nh^{d+2}}\rightarrow 0,\,h\rightarrow 0$. Let $D_{h}$ and $\widehat{D}_{h}$ and $\widehat{D}^{*}_{n}$ be the density level sets with level $\lambda$ for $p_{h}$ and $\widehat{p}_{h}$ and $\widehat{p}_{h}^{*}$. Then

The proof of this Theorem is similar to the proof of Theorem 28, so we omit the details. The basic idea is as follows. By equation (30), the quantile of $\widehat{M}_{n}$ can be used to construct confidence sets. We then show that $\sqrt{nh^{d}}{M}_{n}$ is approximated by the maximum of a Gaussian process (similar to Theorem 3 and made explicit in Chernozhukov et al. 2014a). Finally, we show that the bootstrap $\sqrt{nh^{d}}M^{*}_{n}$ converges to $\sqrt{nh^{d}}{M}_{n}$ as in Theorem 28.

This method embodied in Theorem 5 is very similar to the methods in Jankowski and Stanberry (2012) and Mammen and Polonik (2013). Jankowski and Stanberry (2012) proposes to construct a confidence set of the form

where $A\triangle B=\{x:x\in A,x\notin B\}\cup\{x:x\in B,x\notin A\}$ is the symmetric difference between sets. Then they use the upper quantile of $R_{n}$ to construct a confidence set of a similar form to (29) and apply the bootstrap to estimate the quantile. Their bootstrap consistency relies on Neumann’s method (Proposition 3.1 in Neumann (1998)) and they assume that $h$ converges fast enough so that one can ignore the bias for estimating the original density $p$. Actually, under their assumptions, our proposed bootstrap confidence sets (from both methods 1 and 2) are also consistent for the original level set $D$ since the bias converges faster than the stochastic variation. The method in Mammen and Polonik (2013) should have higher power than our method 2 since they consider taking the supremum over a smaller region.

We may use a variance stabilizing transform to obtain an adaptive confidence set using similar idea to Chernozhukov et al. (2012). The variance of $\widehat{p}_{h}(x)$ is proportional to $p(x)$. Thus, we may use

and set $\widehat{v}_{1-\alpha}\times\sqrt{\widehat{p}_{h}(x)}$ as an adaptive threshold for constructing the confidence set. Namely, the adaptive confidence set is given by

Using the same approach as in the proof to Theorem 5, we can show that $\widehat{C}^{*}_{n,1-\alpha}$ has asymptotically $1-\alpha$ coverage.

The rate $O\left(\left(\frac{\log^{7}n}{nh^{d}}\right)^{1/8}\right)$ may not be optimal. In Chernozukov et al. (2014), they apply a induction technique that gives a rate of order $n^{-1/6}$ for the Gaussian approximation. Despite not being mentioned explicitly in that paper, we believe that similar technique applies to the empirical process. If this is true, the rate in Theorem 3, 28 and 5 can be further refined to $O\left(\left(\frac{\log^{7}n}{nh^{d}}\right)^{1/6}\right)$.

3 Comparing Methods 1 and 2

Despite the fact that method 2 yields a much larger confidence sets, it has some nice properties. First, method 2 is very simple: all we need is to compute the bootstrap distribution of supremum loss. Second, method 2 is connected to the methods proposed in Mammen and Polonik (2013) and Jankowski and Stanberry (2012). Third, the confidence sets produced in method 2 are related to level sets with level $\lambda\pm\widehat{m}_{1-\alpha}$. The last property makes it easy to visualize the confidence sets; see Section 5.3.

Here, we consider two simulated datasets to compare the coverage for confidence sets constructed using Hausdorff loss (method 1), $L_{\infty}$ loss (supremum loss; method 2), and scaled $L_{\infty}$ loss (remark 3).

The first dataset is a three-Gaussian mixture. The data is generated from the following distribution:

where $\phi_{2}(x;\mu,\Sigma)$ is the density to bivariate Gaussian with mean vector $\mu$ and Covariance matrix $\Sigma$. $\mathbf{I}_{2}$ is the $2\times 2$ identity matrix and $\mu_{1}=(0,0)^{T}$, $\mu_{2}=(1,0)^{T}$, and $\mu_{3}=(1.5,0.5)^{T}$. We use density level $\lambda=0.3$ and smoothing parameter $h=0.2$. The corresponding level set $D_{h}$ is the red curve in the left panel of Figure 7. We consider three different sample sizes, $N=500,1000,2500$ and compare the coverage of confidence sets using the three methods. The corresponding coverages are given in Table 1. As can be seen from Table 1, all the three methods have the desire nominal coverage.

The second dataset is a four-mixture dataset from Cadre et al. (2009). The data is generated from the following distribution:

with $\pi_{1}=\pi_{2}=\pi_{3}=\pi_{4}=1/5$ and $\pi_{5}=\pi_{6}=1/10$, and $\mu_{1}=(-0.3,-0.3)^{T},\mu_{2}=(3.0,3.0)^{T},\mu_{3}=\mu_{4}=(0,3)^{T},\mu_{5}=\mu_{6}=(3,0)$, and

We use density level $\lambda=0.05$ and smoothing parameter $h=0.2$. The corresponding level set $D_{h}$ is the red curves in the right panel of Figure 7. Again, we consider three different sample sizes $N=500,1000,2500$ and compare the coverage of confidence sets using the three methods. The corresponding coverages are given in Table 2. It is clear that all the three methods have the desire nominal coverage. Moreover, it can be seen from both Table 1 and 2 that the supremum loss method over covers.

4 Pointwise Hypothesis Tests for Level Sets

The confidence sets developed in the previous section are related to two types of local hypothesis tests. For fixed density level $\lambda$ and an arbitrary point $x$, consider the tests of whether $p(x)$ is greater or less than $\lambda$:

When we only want to test just a few points, we can do local tests for each point and control the family-wise error rate to control the type 1 error rate. Usually, however, we are interested conducting the local test at many or even an infinite number of points (like a region), making it difficult to control type 1 error simultaneously.

Inverting the confidence sets of the previous subsection gives a solution to this problem. Let $\widehat{S}_{n,1-\alpha}$ be a confidence set for $D_{h}$ and let $L_{h}$ and $V_{h}$ be, respectively, the upper and lower lambda level sets. Then the decision rules

In addition, inverting the regions where we cannot reject $H_{\sf in,0}(x)$ and $H_{\sf out,0}(x)$ in (52) yields confidence sets for $L_{h}$ and $V_{h}$. In Figure 6, a $90\%$ confidence regions for the upper level set $L_{h}$ is the union of yellow and blue regions ($T_{\sf out,n}(x)=0$). And a $90\%$ confidence regions for $V_{h}$, the lower level set, is the union of green and blue regions ($T_{\sf in,n}(x)=0$). Thus, with $90\%$ confidence, all yellow regions are above $\lambda$ and the true high density regions should be contained by the yellow and blue regions.

The two local tests described in (50) and (51) are relevant to the problems of level-set clustering (Hartigan, 1975; Polonik, 1995; Rinaldo and Wasserman, 2010; Rinaldo et al., 2012) and anomaly detection (Desforges et al., 1998; Breunig et al., 2000; He et al., 2003; Chandola et al., 2009). Rejecting $H_{\sf in,0}(x)$ can be viewed as evidence that $x$ belongs to a level-set cluster. And a point $x$ where $H_{\sf out,0}(x)$ is rejected can be viewed as anomalous.

Note that one can modify the local testing procedure to control the False Discovery Rate (Benjamini and Hochberg, 1995) rather than familywise error.

Visualization for Multivariate Level Sets

Level-set estimators can reveal useful information about a distribution, but beyond three dimensions, we cannot directly visualize the level sets, making the results difficult to use.

In this section, we propose a novel visualization technique density upper level sets in multidimensions. Any visualization entails some loss of information, but our goals are to preserve important geometric information about the sets, make the overall visualization easy to interpret, and give a method that is efficient to compute. Our method exploits the relationship between level-set clustering and mode clustering (see Figure 9).

A current and commonly used visualization method for level sets is the density tree (Stuetzle, 2003; Klemelä, 2004, 2006; Kent et al., 2013; Balakrishnan et al., 2013). This method considers several density levels, $\lambda_{1}<\cdots<\lambda_{K}$, and computes the number of connected components for the upper density level set at each level. As we increase the density level, some connected components may vanish and others may split into additional components. In the typical case when the underlying density function is a Morse function (i.e., Hessian at critical points is non-degenerate) (Morse, 1925, 1930; Milnor, 1963), the connected component disappears when the density level is above the maximum density value over the component and splits only if the density level passes the density value of some saddle points within the component (Klemelä, 2009). The vanishing and splitting of components as the level changes produces a tree structure. The density tree uses this tree structure as a visualization of the level sets. We refer to Klemelä (2009) for more details.

Density trees display primarily topological information about the underlying density function (Stuetzle, 2003; Kent et al., 2013) but need not preserve or impart other features that may be of interest. Extensions have been proposed that would endow the tree with additional information about the distribution (Klemelä, 2004, 2006, 2009), but this is difficult to do with multiple features and can make the visualization difficult to interpret. See Figure 8 for an example.

In this section, we propose a novel technique that visualizes several density (upper) level sets that preserve some geometric information and that are very easy to understand. Our method is based on the relationship between level set clustering and mode clustering (see Figure 9). Note that in this section, we will focus on density upper level sets.

Our method complements existing tree-based methods. It combine two clustering techniques – level-set and mode clustering – to produce a simple and intuitive visualization.

That is, $\pi_{x}(t)$ starts at $x$ and moves along the gradient of $p$. We define the destination for $\pi_{x}(t)$ as ${\sf dest}(x)=\lim_{t\rightarrow\infty}\pi_{x}(t)$. Let ${\cal M}$ be the collection of all local modes of $p$. It can be shown that ${\sf dest}(x)\in{\cal M}$ except for a set of $x$’s in a set ${\cal B}$ with Lebesgue measure (this set corresponds to the boundaries of clusters). For each mode $m_{j}\in{\cal M}$, we define its basin of attraction as

The regions ${\cal A}_{1},\cdots,{\cal A}_{k}$ are the clusters generated by mode clustering.

Now we recall three facts about an upper density level set $L=\{x:\ p(x)\geq\lambda\}$ (c.f. Figure 9 left and middle panels):

Thus, the upper level sets are covered by the basins of attraction of the local modes (the uncovered regions have Lebesgue measure so we ignore them for visualization).

2 Visualization Algorithm

Given several density levels, $\lambda_{1}<\cdots<\lambda_{K}$, we can overlay the visualization from the previous paragraph from $\lambda_{1}$ to $\lambda_{K}$ to create a “tomographic” visualization of the clusters. This gives a visualization for the density level sets. Figure 10 shows an example for visualizing level sets for a $6$-dimensional and a $10$-dimensional simulation datasets at different density levels. This dataset is from Chen et al. (2014c).

3 Visualizing Level Sets with Confidence

A modified version of our algorithm allows us to visualize multivariate confidence sets for an upper level set at a given level $\lambda$. In particular, we visualize the confidence set for the level sets produced by Method 2 (supremum loss, see Section 4.2) . The advantage of Method 2 here is that this confidence set under different $\alpha$’s is just the density level set at different $\lambda$’s. This makes it easy to visualize $K$ distinct confidence sets for pre-specified levels $\alpha_{1}<\cdots<\alpha_{K}$.

and use the visualization algorithm to create a tomographic visualization. Figure 11 provides an example for visualizing confidence sets for upper level set using the 6 dimensional simulation data in Figure 10.

At the cost of additional computation, we can also use Method 1 (Hausdorff loss) instead of Method 2, combining it with the basins of attractions to visualize the confidence sets. We use method 1 to construct the inner/outer confidence sets and then we find the connected components and partition it by the basins of attraction for local modes and use multidimensional scaling to visualize it in low dimensions.

Discussion

In this paper, we derived the limiting distribution for smoothed density level sets under Hausdorff loss. This result immediately allows us to construct confidence sets for the smoothed level set. We developed two bootstrapping methods to construct the confidence sets, and we showed that both methods are consistent. These confidence sets can be inverted to construct multiple local tests of whether a point’s density value is above or below a given level, which has application to related problems such as anomaly detection. Finally, we developed a new visualization method that is informative and interpretable, even in multidimensions.

Although we focused on density level sets in this paper, our methods, including confidence sets and visualization, can be applied to a more general class of problems such as kernel classifier, two-sample tests, and generalized level sets (See Mammen and Polonik 2013 for more details.).

References

Proofs

Assume (K1–2), then for each $t>0$ there exists some $n_{0}$ such that whenever $n>n_{0}$, we have

Proof for Lemma 1. We first prove the lower bound for $\operatorname{{\sf reach}}(D_{h})$ and then we will prove the additional assertions.

Part 1: Lower bound on reach. We prove this by contradiction. Take $x$ near $D_{h}$ such that

We assume that $x$ has two projections onto $D_{h}$, denoted as $b$ and $c$.

Since $b,c\in D_{h}$, $p_{h}(b)-\lambda=p_{h}(c)-\lambda=0$ so that $p_{h}(b)-p_{h}(c)=0$. Now by Taylor’s theorem

Since both $b$ and $c$ are projection points from $x$ onto $D_{h}$,

Recall that $d(x,D_{h})\leq\frac{g_{0}}{\|p_{h}\|_{2,\max}}$ and by Taylor’s theorem,

so that $|t_{b}|\|p_{h}\|_{2,\max}<1$. Note that the lower bound $g_{0}$ in the last inequality is because $d(x,D_{h})<\frac{\delta_{0}}{2}$ so it follows from assumption (G). Plugging in this result into the last equality of (60), we conclude that $\|b-c\|=0$. This shows $b=c$ so that we have a unique projection. Thus, whenever $d(x,D_{h})<\left(\frac{\delta_{0}}{2},\frac{g_{0}}{\|p_{h}\|^{*}_{2,\max}}\right)$, we have a unique projection onto $D_{h}$ and thus we have proved the lower bound on reach.

Part 2: The three assertions. The first assertion is trivially true when $\|p_{h}-q\|^{*}_{2,\max}$ is sufficiently small since assumption (G) only involves gradients (first derivatives).

The second assertion follows from the lower bound on reach. By assertion 1, (G) holds for $q$. And the lower bound on reach is bounded by gradient and second derivatives so that we have the prescribed bound.

The third assertion follows from Theorem 1 in Chazal et al. (2007) which states that if two $d-1$ dimensional smooth manifolds $M_{1}$ and $M_{2}$ have Hausdorff distance being less than $(2-\sqrt{2})\min\{\operatorname{{\sf reach}}(M_{1}),\operatorname{{\sf reach}}(M_{2})\}$, then $M_{1}$ and $M_{2}$ are normal compatible to each other. Now by Theorem 8, the Hausdorff distance between $D_{h}$ and $D(q)$ is at rate $O(\|p_{h}-q\|_{1,\max})$ so that this assertion is true when $\|p_{h}-q\|_{2,\max}$ is sufficiently small. $\square$

Proof of Lemma 2. Let $x\in D_{h}$. We define $\Pi(x)\in D_{h}$ to be the projected point onto $\widehat{D}_{h}$. By Lemma 1 and Theorem 8, when $\|\widehat{p}_{h}-p_{h}\|^{*}_{2,\max}\rightarrow 0$, ${\sf Haus}(D_{h},\widehat{D}_{h})\overset{P}{\rightarrow}0$ so that $\Pi(x)$ is unique. Thus, we assume $\Pi(x)$ is unique.

Now since $\Pi(x)\in\widehat{D}_{h}$ and $x\in D_{h}$, $\widehat{p}_{h}(\Pi(x))-p_{h}(x)=0$. Thus, by Taylor’s theorem

Note that $x-\Pi(x)$ is normal to $\widehat{D}_{h}$ at $\Pi(x)$ so that it points toward the same direction as $\nabla\widehat{p}_{h}(\Pi(x))$. Thus, (62) can be rewritten as

By Taylor’s theorem, $\nabla\widehat{p}_{h}(\Pi(x))$ is close to $\nabla p_{h}(x)$ in the sense that

In addition, $O(\|x-\Pi(x)\|))$ is bounded by $O({\sf Haus}(\widehat{D}_{h},D_{h}))$ which is at rate $O(\|\widehat{p}_{h}-p_{h}\|^{*}_{1,\max})$ due to Theorem 8. Putting this together with (63), we conclude

Note that the left hand side can be written as

This holds uniformly for all $x\in D_{h}$ and note that the definition of ${\cal F}$ is

Proof for Theorem 3. The proof for Theorem 3 follows the same procedure as the proof of Theorem 6 in Chen et al. (2014b). The proof contains two parts: Gaussian approximation and anti-concentration.

Part 1: Gaussian approximation. Basically, we will show that

First, when $\|\widehat{p}_{h}-p_{h}\|$ is sufficiently small, $\widehat{D}_{h}$ and $D_{h}$ are normal compatible to each other by Lemma 1. Then by the property of normal compatible,

Note that this result basically follows from the same derivation of Proposition 3.1 in Chernozhukov et al. (2014c) with the fact that $g\equiv 1$ in their definition.

Combining equations (70) and (71) and pick $t=1/\sqrt{nh^{d+2}}$, we have that for $n$ is sufficiently large and $\gamma\in(0,1)$,

Part 2: Anti-concentration. To obtain the desired Berry-Esseen bound, we apply the anti-concentration inequality in Chernozhukov et al. (2014c) and Chernozhukov et al. (2014a).

From Lemma 10 and equation (72), there exists some constant $A_{6}$ such that

Now pick $\gamma=\left(\frac{\log^{7}n}{nh^{d}}\right)^{1/8}$ and use the fact that $\frac{1}{\sqrt{nh^{d+2}}}$ and $2e^{-\sqrt{nh^{d+2}}A_{2}}$ converges faster than the other terms; we obtain the desired rate. $\square$

Proof for Theorem 4. This proof follows the same strategy for the proof of Theorem 7 in Chen et al. (2014b). We prove the Berry-Esseen type bound first and then show that the coverage is consistent. We prove the Berry-Esseen bound in two simple steps: Gaussian approximation and support approximation.

Thus, if we sample from $\widehat{p}_{h}$ and consider estimating $\widehat{p}_{h}$ by $\widehat{p}^{*}_{h}$, we are doing exactly the same procedure of estimating $p_{h}$ by $\widehat{p}_{h}$. Therefore, Lemma 2 and Theorem 3 hold for approximating ${\sf Haus}(\widehat{D}^{*}_{h},\widehat{D}_{h})$ by a maxima for a Gaussian process. The difference is that the Gaussian process is defined on

since the “parameter (level sets)” being estimated is $\widehat{D}_{h}$ (the estimator is $\widehat{D}^{*}_{h}$). Note that ${\cal F}_{n}$ is very similar to ${\cal F}$ except the denominator is slightly different and the support $\widehat{D}_{h}$ is also different from $D_{h}$. That is, we have

Step 2: Support approximation. In this step, we will show that

The first approximation can be shown by using the Gaussian comparison lemma (Theorem 2 in Chernozhukov et al. (2014b); also see Lemma 17 in Chen et al. (2014b)). We do the same thing as Step 3 in the proof of Theorem 8 in Chen et al. (2014b) so we omit the details. Essentially, given any $\epsilon>0$, we can construct a pair of balanced $\epsilon$-nets for both ${\cal F}$ and ${\cal F}_{n}$, denoted as $\{g_{1},\cdots,g_{K}\}$ and $\{g^{n}_{1},\cdots,g^{n}_{K}\}$ so that $\max_{j}\|g_{j}-g^{n}_{j}\|^{*}_{\max}=O(\|\widehat{p}_{h}-p_{h}\|^{*}_{1,\max})$. Then this $\epsilon$-net leads to

Now comparing the above result to Theorem 3 and using the fact that the first big-O term dominates the second term (the first is of rate $-1/8$ for $n$ but the second term is at rate $-1/6$ by Theorem 9), we conclude the result for first assertion.

For the coverage, let $W_{n}={\sf Haus}(\widehat{D}_{h},D_{h})$ and $w_{n,1-\alpha}=F^{-1}_{W_{n}}(1-\alpha)$. Since $D_{h}\subset\widehat{D}_{h}\oplus{\sf Haus}(\widehat{D}_{h},D_{h})$, we have

Now by the first assertion, the difference for $w_{n,1-\alpha}$ and the bootstrap estimate $w^{*}_{n,1-\alpha}$ differs at rate $O\left(\left(\frac{\log^{7}n}{nh^{d}}\right)^{1/8}\right)$, which completes the proof.