Nonparametric Confidence Regions for Level Sets: Statistical Properties and Geometry

TL;DR

This paper explores the theoretical properties and geometric aspects of constructing nonparametric confidence regions for density level sets, supported by numerical studies.

Contribution

It provides a comprehensive analysis of both vertical and horizontal variation methods for confidence regions, including their geometric relationships and finite sample performance.

Findings

Theoretical insights into the behavior of confidence regions via large sample theory

Analysis of geometric and topological influences on finite sample performance

Numerical studies validating the theoretical and geometric discussions

Abstract

This paper studies and critically discusses the construction of nonparametric confidence regions for density level sets. Methodologies based on both vertical variation and horizontal variation are considered. The investigations provide theoretical insight into the behavior of these confidence regions via large sample theory. We also discuss the geometric relationships underlying the construction of horizontal and vertical methods, and how finite sample performance of these confidence regions is influenced by geometric or topological aspects. These discussions are supported by numerical studies.

Figures

Tables

Table 1. Table 1: Simulation results

		Case 1			Case 2			Case 3			Case 4
		C.P.	$\overline{\lambda }$	$\overline{P}$	C.P.	$\overline{\lambda }$	$\overline{P}$	C.P.	$\overline{\lambda }$	$\overline{P}$	C.P.	$\overline{\lambda }$	$\overline{P}$
n=200	H	0.98	8.40	0.63	1.00	13.18	0.82	1.00	41.21	0.49	1.00	59.46	0.88
	V.e	0.87	6.91	0.53	0.91	6.87	0.53	0.95	35.16	0.15	0.90	12.17	0.47
	V	0.91	7.33	0.56	0.94	7.28	0.56	0.95	36.17	0.17	0.85	13.14	0.49
	V.bc	0.50	6.28	0.46	0.56	6.26	0.46	0.66	20.63	0.15	0.53	11.75	0.43
	V.us	0.94	13.55	0.77	0.96	13.45	0.77	0.96	28.58	0.26	0.88	29.82	0.71
n=1000	H	0.95	4.26	0.34	0.98	7.75	0.6	1.00	21.22	0.18	0.98	33.35	0.76
	V.e	0.90	3.77	0.30	0.89	3.79	0.30	0.92	36.31	0.09	0.9	7.41	0.31
	V	0.94	4.09	0.33	0.94	4.11	0.33	0.92	37.37	0.11	0.89	8.00	0.33
	V.bc	0.61	3.54	0.27	0.61	3.55	0.28	0.64	25.60	0.1	0.63	7.41	0.31
	V.us	0.94	6.27	0.47	0.93	6.30	0.47	0.96	30.81	0.15	0.94	13.78	0.49
n=5000	H	0.93	2.32	0.19	0.97	4.47	0.36	0.98	9.67	0.07	0.99	15.29	0.58
	V.e	0.90	2.17	0.17	0.88	2.16	0.17	0.94	9.38	0.06	0.90	4.77	0.22
	V	0.94	2.38	0.19	0.96	2.37	0.19	0.96	13.51	0.07	0.92	5.17	0.23
	V.bc	0.74	2.08	0.16	0.73	2.07	0.16	0.74	8.82	0.06	0.71	4.84	0.22
	V.us	0.97	3.52	0.28	0.96	3.49	0.27	0.97	33.49	0.10	0.95	7.89	0.33
	V.ls	1.00	5.96	0.44	1.00	6.39	0.47	1.00	34.45	0.13	1.00	13.85	0.50
n=50000	H	0.93	1.04	0.08	0.92	2.06	0.17	0.94	3.86	0.03	1.00	7.89	0.36
	V.e	0.92	1.01	0.08	0.92	0.98	0.08	0.90	3.68	0.03	0.93	2.73	0.13
	V	0.96	1.11	0.09	0.95	1.09	0.09	0.94	4.65	0.03	0.97	2.95	0.14
	V.bc	0.81	0.98	0.08	0.84	0.96	0.08	0.77	3.57	0.03	0.83	2.75	0.13
	V.us	0.97	1.63	0.13	0.96	1.60	0.13	0.96	6.62	0.05	0.98	4.09	0.19
	V.ls	1.00	1.64	0.13	1.00	1.72	0.14	1.00	6.44	0.05	1.00	4.09	0.19

Equations

$$\mathcal{L}=\{x\in {\mathbb{R}}^d:f(x)\ge c\}\text{and}\mathcal{M}=\{x\in {\mathbb{R}}^d:f(x)=c\},$$

$$\displaystyle\widehat{f}(x)=\frac{1}{nh^{d}}\sum_{i=1}^{n}K\Big{(}\frac{X_{i}-x}{h}\Big{)},\;\;x\in\mathbb{R}^{d},$$

$$\hat{\mathcal{L}}=\{x\in {\mathbb{R}}^d:\hat{f}(x)\ge c\}\text{and}\hat{\mathcal{M}}=\{x\in {\mathbb{R}}^d:\hat{f}(x)=c\}.$$

$$\hat{\beta }(x)=\frac{1}{2}{h}^2\int u_1^{2}K(u)du\sum _{j=1}^d\frac{{\partial }^2}{\partial x_j\partial x_j}{\hat{f}}_l(x),$$

$${\mathcal{L}}^E=\{x\in {\mathbb{R}}^d:\mathbb{E}\hat{f}(x)\ge c\}\text{and}{\mathcal{M}}^E=\{x\in {\mathbb{R}}^d:\mathbb{E}\hat{f}(x)=c\},$$

$$\displaystyle\widehat{C}=\widehat{f}^{-1}\big{[}c-\widehat{a}_{n},c+\widehat{a}_{n}\big{]}.$$

$$\hat{C}=\{y\in {\mathbb{R}}^d:\Vert y-x\Vert \le {\hat{b}}_n(x),x\in {\hat{f}}^{-1}(c)\},$$

$${\beta }_{n,\gamma }^{(k)}={\gamma }^2+\sqrt{\frac{\log n}{n{\gamma }^{d+2k}}},\text{and}{\beta }_{n,\gamma }^{(k),E}=\sqrt{\frac{\log n}{n{\gamma }^{d+2k}}},k=0,1,2,3.$$

$$\hat{b}(\alpha )$$

$$\displaystyle\hskip 113.81102pt+\log\bigg{\{}\frac{(2d-2)^{d/2-1}s_{K}^{d-1}}{\sqrt{2}\pi^{d/2}}\widehat{\mathscr{V}_{d-1}(\mathcal{M})}\bigg{\}}\bigg{]}.$$

$$\displaystyle\widehat{a}^{(d)}_{1-\alpha}=\begin{cases}\frac{\widehat{b}(\alpha)\sqrt{\|K\|_{2}^{2}c}}{\sqrt{nh^{d}}}&\text{for }\;d\geq 2\\ \frac{\Phi^{-1}\big{(}(1-\alpha)^{1/\widehat{N}}\big{)}\sqrt{\|K\|_{2}^{2}c}}{\sqrt{nh^{d}}}&\text{for }\;d=1,\end{cases}$$

$${\hat{C}}_{n,1}(1-\alpha )=({\hat{f}}^{bc}{)}^{-1}\left[c-{\hat{a}}_{1-\alpha }^{(d)},c+{\hat{a}}_{1-\alpha }^{(d)}\right].$$

$$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\mathcal{M}\subset\widehat{C}_{n,1}(1-\alpha)\Big{)}=1-\alpha.$$

$$\displaystyle\widehat{C}_{n,1}^{-}(1-\alpha)=(\widehat{f}^{\,bc})^{-1}\Big{[}c-\widehat{a}^{(d)}_{1-\alpha},+\infty\Big{)}$$

$$\displaystyle\widehat{C}_{n,1}^{+}(1-\alpha)=(\widehat{f}^{\,bc})^{-1}\Big{[}c+\widehat{a}^{(d)}_{1-\alpha},+\infty\Big{)}.$$

$$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\widehat{C}_{n,1}^{+}(1-\alpha)\subset\mathcal{L}\subset\widehat{C}_{n,1}^{-}(1-\alpha)\Big{)}=1-\alpha.$$

$$\hat{\mathcal{L}}={\mathrm{argmin}}_{A\in \mathcal{A}}[P_n(A)-c\lambda (A)].$$

$$T(D_{n})=\sup_{x\in D_{n}}\big{|}\sqrt{nh^{d}\,}\,(\widehat{f}(x)-f(x))\big{|},\text{ \rm or }T^{E}(D_{n})=\sup_{x\in D_{n}}\big{|}\sqrt{nh^{d}\,}\,(\widehat{f}(x)-{\mathbb{E}}\widehat{f}(x))\big{|},$$

$${\hat{C}}_{n,2}^{\ast }(1-\alpha )={\hat{f}}^{-1}[c-{\hat{c}}_{1-\alpha }^{\ast },c+{\hat{c}}_{1-\alpha }^{\ast }]$$

$${\hat{C}}_{n,2}^{\ast ,E}(1-\alpha )={\hat{f}}^{-1}[c-{\hat{c}}_{1-\alpha }^{\ast ,E},c+{\hat{c}}_{1-\alpha }^{\ast ,E}].$$

$${\hat{C}}_{n,2}^{\ast ,-}(1-\alpha )={\hat{f}}^{-1}[c-{\hat{c}}_{1-\alpha }^{\ast },+\infty ),{\hat{C}}_{n,2}^{\ast ,+}(1-\alpha )={\hat{f}}^{-1}[c+{\hat{c}}_{1-\alpha }^{\ast },+\infty ),$$

$${\hat{C}}_{n,2}^{\ast ,E,-}(1-\alpha )={\hat{f}}^{-1}[c-{\hat{c}}_{1-\alpha }^{\ast ,E},\infty ),\text{and}{C}_{n,2}^{\ast ,E,+}(1-\alpha )={\hat{f}}^{-1}[c+{\hat{c}}_{1-\alpha }^{\ast ,E},+\infty ).$$

$${\Psi }_n(\gamma )=\gamma \sqrt{n{h}^d\log n}+h\log n+\log n\left(\sqrt{{\beta }_{n,h}^{(0),E}}+\sqrt{h^d\log n}\right),\gamma >0.$$

$${\gamma }_n^E$$

$${\gamma }_n$$

$$\displaystyle\mathbb{P}\left(\mathcal{M}^{E}\subset\widehat{C}_{n,2}^{*,E}(1-\alpha)\right)=(1-\alpha)+O\big{(}\Psi_{n}(\gamma_{n}^{E})\big{)},\quad\text{and}$$

$$\displaystyle\mathbb{P}\left(\widehat{C}_{n,2}^{*,E,+}(1-\alpha)\subset\mathcal{L}^{E}\subset\widehat{C}_{n,2}^{*,E,-}(1-\alpha)\right)=(1-\alpha)+O\big{(}\Psi_{n}(\gamma_{n}^{E})\big{)}.$$

$$\displaystyle\mathbb{P}\left(\mathcal{M}\subset\widehat{C}_{n,2}^{*}(1-\alpha)\right)=(1-\alpha)+O\big{(}\Psi_{n}(\gamma_{n})\big{)},\quad\text{and}$$

$$\displaystyle\mathbb{P}\left(\widehat{C}_{n,2}^{*,+}(1-\alpha)\subset\mathcal{L}\subset\widehat{C}_{n,2}^{*,-}(1-\alpha)\right)=(1-\alpha)+O\big{(}\Psi_{n}(\gamma_{n})\big{)}.$$

$${\rho }_n^E=O(h(\log n{)}^{3/2})=O(n^{-1/(d+4)}(\log n{)}^{3/2}).$$

$${\rho }_n=O\left(n^{-1/(d+4)}(\log n{)}^{5/4}+{\beta }_{n,g}^{(2)}\sqrt{\log n}+\sqrt{{\beta }_{n,g}^{(1)}}\log n\right).$$

$${\hat{C}}_{n,3}^{\ast }(1-\alpha )=({\hat{f}}^{bc}{)}^{-1}[c-{\hat{c}}_{1-\alpha }^{\ast ,E},c+{\hat{c}}_{1-\alpha }^{\ast ,E}].$$

$${\hat{C}}_{n,3}^{\ast ,-}(1-\alpha )=({\hat{f}}^{bc}{)}^{-1}[c-{\hat{c}}_{1-\alpha }^{\ast ,E},+\infty ),\text{ and }$$

Full text

62G20,

62G05,

Extreme value distribution,

level sets,

nonparametric surface estimation,

integral curves,

kernel density estimation,

\startlocaldefs\endlocaldefs

t2Partially supported by NSF grant DMS 1821154. t3Partially supported by NSF grant DMS 1107206.

1. 1 Introduction

2. 2 Confidence regions based on vertical variation

3. 2.1 Confidence regions based on asymptotic distribution

4. 2.2 Bootstrap confidence regions

5. 2.2.1 A bootstrap confidence region based on explicit bias correction

6. 3 Confidence regions based on horizontal variation

7. 3.1 Methods based on integral curves

8. 3.1.1 Confidence regions for ${\cal M}$ and ${\cal L}$ using local adjustment by the gradient

9. 3.1.2 Confidence regions for ${\cal M}^{E}$ and ${\cal L}^{E}$ adjusted by gradient

10. 3.1.3 Confidence regions not locally adjusted by the gradient

11. 3.2 Horizontal variation based methods not based on integral curves

12. 4 Performance of confidence regions and geometry

13. 5 Simulations

14. 6 Conclusion

15. 7 Proofs

16. 7.1 Proof of Theorem 2.1

17. 7.2 Proof of Corollary 2.1

18. 7.3 Proof of Theorem 2.2

19. 7.4 Proof of Theorem 2.3

20. 7.5 Proofs for Section 3

21. 8 Appendix

1 Introduction

For a density function $f$ on $\mathbb{R}^{d},\;d\geq 1,$ we define the superlevel set of $f$ at level $c$ and the corresponding isosurface or contour as

[TABLE]

respectively. The dependence on the level $c$ is suppressed in our notation, for $c$ is fixed throughout the manuscript. We will study methods for constructing confidence regions for ${\cal M}$ and ${\cal L}$ based on an iid sample $X_{1},X_{2},\cdots,X_{n}$ from $f$. While new methods for constructing confidence regions are proposed below, and a comparison of these and other existing methods is provided, the overarching goal of this work is to provide a critical discussion of the advantages and disadvantages of the various existing methods. This will, among others, also lead to insight about the influence of geometry on the statistical performance of the confidence regions.

The estimation of level sets (or isosurfaces) has received quite some interest in the literature. For some earlier work in the context of density level set estimation see Hartigan (1987), Polonik (1995), Cavalier (1997), Tsybakov (1997), Walther (1997). There are relations to density support estimation, (e.g. Cuevas et al. 2004, Cuevas 2009), clustering (e.g. Cuevas et al. 2000, Rinaldo et al. 2010), classification (e.g., Mammen and Tsybakov, 1999, Hall and Kang, 2005, Steinwart et al. 2005, Audibert and Tsybakov 2007), anomaly detection (e.g. Breuning et al. 2000, Hodge and Austin 2004), and more. Bandwidth selection for level set estimation is considered by Samworth and Wand (2010), and Qiao (2018a). Applications of level set estimation exist in many fields, such as astronomy (Jang 2006), medical imaging (Willett and Nowak, 2005), geoscience (Sommerfeld et al. 2015), and others. See Mason and Polonik (2009) and Mammen and Polonik (2013) for further literature. Level sets of functions also play a central role in topological data analysis, in particular in ‘persistent homology’, where the topological properties of level sets as a function of the parameter $c$ are used for statistical analysis; e.g. see Fasy et al. (2014).

All the methods discussed in this paper are based on a kernel estimator

[TABLE]

and its derivatives. Here $K$ is a $d$-dimensional kernel, and $h>0$ is the bandwidth. As estimators for the superlevel sets or the corresponding contours, we are considering plug-in estimators given by

[TABLE]

It is well-known that $\widehat{f}(x)$ is a biased estimator. For a twice continuously differentiable kernel $K$, we will also consider a de-biased version of the kernel estimator, given by $\widehat{f}^{\,bc}=\widehat{f}-\widehat{\beta}$ with

[TABLE]

where $\widehat{f}_{l}$ denotes a kernel density estimator using bandwidth $l$, which can be different from $h$. See Chen (2017) for a recent study of the de-biased estimator.

In some of the recent related literature (see, e.g., Chen et al., 2017), the bias is ignored entirely, and the target is redefined as a ‘smoothed’ version of the superlevel set, or its contour, given by

[TABLE]

respectively, where here and in what follows, we use the superscript $E$ to indicate that the definition is based on an ‘expected quantity’. We will also consider this target for some of our methods.

For an interval $[a,b]\subset{\mathbb{R}},$ and a real valued function $g$, we denote by $g^{-1}[a,b]$ the pre-image of $[a,b]$ under $g$. For $a=b$, we also use the standard notation $g^{-1}(a)$ for the pre-image. Two different types of confidence regions are considered. One of them is based on vertical variation. Such confidence regions for ${\cal M}$ or ${\cal M}^{E}$ are of the form

[TABLE]

Corresponding lower and upper confidence regions for ${\cal L}$ (and ${\cal L}^{E}$) will be of the form \widehat{f}^{-1}\big{[}c+\widehat{a}_{n},\infty\big{)} and \widehat{f}^{-1}\big{[}c-\widehat{a}_{n},\infty\big{)}, respectively. The crucial question then is, how to choose the quantity $\widehat{a}_{n}$ in order to achieve a good performance. The other type of confidence region considered here is based on horizontal variation. Such confidence regions for ${\cal M}$ (or ${\cal M}^{E}$) are of the form

[TABLE]

where the random variable $\widehat{b}_{n}(x)$ either does not depend on $x$, i.e. $\widehat{b}_{n}(x)=\widehat{b}_{n}$, or $\widehat{b}_{n}(x)=\widehat{a}_{n}/\|\nabla\widehat{f}(x)\|$, for some random variable $\widehat{a}_{n}$ not depending on $x$. Once horizontal variation based confidence regions $\widehat{C}$ for ${\cal M}$ are constructed, we can obtain corresponding confidence regions for ${\cal L}$: Upper bounds (sets) are obtained by adding $\widehat{C}$ to $\widehat{\cal L}$, and lower bounds are constructed by subtraction. Confidence regions for ${\cal M}^{E}$ and ${\cal L}^{E}$ are constructed similarly.

For a simple heuristic underlying the choice of $\widehat{b}_{n}(x)=\widehat{a}_{n}/\|\nabla\widehat{f}(x)\|$, consider the one-dimensional case. For small $\widehat{a}_{n}$, we have $|\widehat{a}_{n}|\approx|\widehat{f}(x+\widehat{b}_{n}(x))-\widehat{f}(x)|=|\widehat{f}(x+\widehat{b}_{n}(x))-c|$ for $x\in\widehat{f}^{-1}(c)$, so that $\widehat{a}_{n}$ reflects vertical variation, while $\widehat{b}_{n}(x)$ represents horizontal variation. For more details see Section 3.

Geometrically, confidence regions based on horizontal variation as in (1.2), but with a constant $\widehat{b}_{n},$ consist of tubes of constant width put around the estimated contour $\widehat{\cal M}.$ In other words, the maximum (horizontal) distance to ${\cal M},$ or ${\cal M}^{E},$ is being controlled. Confidence regions based on vertical variation as in (1.1), as well as the regions based on horizontal variation with non-constant $\widehat{b}_{n}(x),$ have varying width, which contain information about the slope of the density. Different constructions of $\widehat{a}_{n}$ and $\widehat{b}_{n}$ will be discussed. One of them is based on extreme value theory for Gaussian fields indexed by manifolds, and the others on various bootstrap approximations. Some of the constructions of the horizontal methods also involve the estimation of integral curves.

The same type of confidence regions as discussed above, can also be constructed with $\widehat{f}(x)$ replaced by the de-biased density estimator $\widehat{f}^{bc}(x)$. We note that in the literature, different approaches to the removal of bias effect in confidence bands or regions using kernel-type estimators. One is based on undersmoothing of the original kernel density estimator. This approach is not considered in some detail here (but see Remark 2.2 below). Alternatively, one can use explicit bias corrections, or the smoothed bootstrap. Both of these approaches are being discussed in our work. For further relevant recent literature see Chen (2017) and Calonico et al. (2018a).

Bootstrap confidence regions for a density superlevel set and/or isosurface based on vertical variation can be found in Mammen and Polonik (2013) and Chen et al. (2017). The latter also constructed a confidence set based on horizontal variation, and such constructions can also be found in Chen (2017). The confidence sets proposed here are compared to these methods. In a different setting Jankowski and Stanberry (2014) consider confidence regions for the expected value of a random sets (or its boundary) based on repeated observations of the random set.

In summary, the contributions of the manuscript are to

• (i)

derive asymptotically valid confidence regions based on vertical variation using extreme value distributions of kernel estimators indexed by certain manifolds;

• (ii)

use bootstrap methods to construct confidence regions based on vertical variation in order to improve finite sample performance of the confidence regions in (i);

• (iii)

derive asymptotically valid bootstrap confidence regions based on horizontal variation;

• (iv)

discuss geometric connections between the different constructions of confidence regions;

• (v)

provide numerical studies to compare the finite sample performance of the various confidence regions developed here and in the literature;

• (vi)

critically discuss advantages and disadvantages of the two different types of constructing confidence regions (horizontal and vertical variation).

The theoretical result underlying (i), see Theorem 2.1, provides a closed form of asymptotically valid confidence regions for superlevel sets and isosurfaces. To the best of our knowledge, this is the first result of this type, and it provides important qualitative insight into the underlying problem. Due to the well-known slow convergence properties of extreme value distributions, the construction of bootstrap confidence regions appear of higher practical relevance.

In a recent paper, Qiao and Polonik (2018) derived the asymptotic distribution of extrema of certain rescaled non-homogeneous Gaussian fields. This result provides an important tool for the construction of our confidence regions based on large sample distribution theory. The paper Qiao and Polonik (2016) on ridge (or filament) estimation is also a source of inspiration for the work presented here.

The paper is organized as follows. Section 2 is considering vertical variation based methods. We first present a result on the behavior of the coverage probability of asymptotic confidence regions for isosurfaces and levels sets, and then we construct a bootstrap-based confidence region. Section 3 is dedicated to horizontal methods. Section 3.1 presents the construction of two confidence regions of the form (1.2) using integral curves. Section 3.2 discusses a horizontal bootstrap based confidence region related to the Hausdorff-distance based method of Chen et al. (2017). In Section 4 we discuss how the finite sample behavior of our methods is influenced by some geometric or topological properties of estimated superlevel sets etc. Simulation results presented in Section 5 compare the various methods. Section 6 presents some concluding discussions. The proofs of the technical results are presented in Section 7.

2 Confidence regions based on vertical variation

The following notation is used throughout the manuscript. For a sequence $\gamma>0,$ we let

[TABLE]

For a kernel density estimator based on a bandwidth $h,$ the quantity $\beta_{n,h}^{(0)}$ equals the rate of the uniform deviation from the density under standard assumptions (satisfied in our setting). In other words, we have uniform consistency of the kernel estimator if $\beta_{n,h}^{(0)}\to 0$ as $n\to\infty.$ The quantities $\beta_{n,h}^{(k)},k=1,2,3$ have the same interpretation when considering the kernel estimator of the $k$-th derivatives of the density. Similarly, $\beta_{n,h}^{(k),E}$ is the standard uniform rate of convergence of the kernel density estimator of the $k$-th derivative when centered at its expectation, $k=0,\ldots,3$.

2.1 Confidence regions based on asymptotic distribution

Our first result provides asymptotically valid confidence regions for the isosurface $\mathcal{M}$ and the superlevel set $\mathcal{L}$. Before formulating the theorem we introduce the underlying assumptions and some more notation.

Assumptions:

(F1) The probability density $f$ has bounded, continuous derivatives up to fourth order. There exists $m>0$ such that $\int\|x\|^{m}f(x)dx<\infty$.

(F2) There exist $\delta_{0}>0$ and $\epsilon_{0}>0$ such that $\|\nabla f(x)\|>\epsilon_{0}$ for $x\in\{x:c-\delta_{0}\leq f(x)\leq c+\delta_{0}\}.$

(K) The kernel function $K:\;\mathbb{R}^{d}\rightarrow\mathbb{R}$ is symmetric about zero, with support contained in $[-1,1]^{d}$, and is continuously differentiable up to fourth order.

(A) Both $f$ and $K$ are continuously differentiable up to $d$-th order.

(H1)k. $h\rightarrow 0$, and $\beta_{n,h}^{(k),E}\rightarrow 0$ as $n\rightarrow\infty$, where $k=0$ or $2$.

(H2)k The bandwidth $l$ used for bias correction, satisfies $l\rightarrow 0$ and $\beta_{n,l}^{(2k),E}\rightarrow 0$, where $k=0$ or $2$. In addition, $(h/l)^{\frac{d+4}{2}}\log{n}\rightarrow 0$ and $\sqrt{nh^{d+4}\log{n}\,}\,l^{2}\rightarrow 0$.

Remark 2.1

a) Assumption (*F1) can be weakened when only considering confidence regions for the smoothed isosurface ${\cal M}^{E}$. Only continuous second order derivatives of $f$ are needed in this case. For simplicity, we just use the stronger assumption (F1) throughout the manuscript.

b) Assumptions (F1) (even the weakened version discussed in 1.) and (F2) together imply that the $(d-1)$-dimensional isosurface $\mathcal{M}$ has positive reach (e.g. see Lemma 1 of Chen et al., 2017), which is a necessary condition for applying the result in Qiao and Polonik (2018). Assumptions on the reach (introduced by Federer 1959) are used in a number of studies of geometric properties of manifolds, etc. Also note that assumption (F2) implies that the level $c>0$.

c) It seems likely that the assumption of a non-negative kernel can be replaced by a higher-order kernel. This is because our focus is on the superlevel sets (or isosurfaces) with $c>0$ and regions there the density estimator is negative is irrelevant to our estimation. Even for the bootstrap method, technically we can bootstrap from a truncated and normalized density estimator. However, the truncation might incur further technical considerations.

d) The smoothness requirements for $f$ and $K$ in assumption (A) are needed to enable to use of the Rosenblatt transform (see Rosenblatt, 1976) in the proof of Theorem 2.1. It will not be needed for bootstrap based methods. This is why we did not combine assumptions (A), (F1) and (K).

e) Assumptions (H1)0 and (H2)0, both very mild, are used in large sample confidence regions based on the vertical method. They guarantee the uniform consistency of the de-biased estimator $\widehat{f}^{bc}=\widehat{f}-\widehat{\beta}$. Specifically, (H1)0 is used for $\widehat{f}$ centered at its expectation, while (H2)0 is used for the bias correction part $\widehat{\beta}$. We use stronger assumptions (H1)2 and (H2)2 for confidence regions based on the horizontal method. In particular, these assumptions guarantee the uniform consistency of the second derivatives of $\widehat{f}^{bc}$.

f) Under our assumptions, in the case of $d=1$, $\cal M$ is a union of finitely many points, say $N$. Denote $\mathcal{M}=\{x_{i},i=1,\cdots,N\}$. **

We need to introduce further notation. Let $\beta(x)=\mathbb{E}\widehat{f}(x)-f(x)$ be the bias. Set $\|K\|_{2}^{2}=\int K^{2}(u)du,$ and let $\mathscr{V}_{d-1}$ denote $(d-1)$-dimensional Hausdorff measure. Let $\widehat{\mathscr{V}_{d-1}(\mathcal{M})}$ be an estimator of $\mathscr{V}_{d-1}(\mathcal{M})$. Some specific estimators will be given in Remark 2.2c). For $d\geq 2$ and s^{2}_{K}=\frac{\int\big{[}\frac{\partial}{\partial u_{1}}K(u)\big{]}^{2}du}{2\int K^{2}(u)du}, we set

[TABLE]

where $0<\alpha<1$, and $z(\alpha)=-\log[-\frac{1}{2}\log(1-\alpha)]$. The quantity $\widehat{b}(\alpha)$ is based on the extreme value behavior of a Gaussian field indexed by ${\cal M}$ (see (7.1) given in Theorem 7.1). We further need the following, which, for $d\geq 2$, simply is a scaled version of $\widehat{b}(\alpha):$

[TABLE]

where $\Phi$ is the standard normal cdf, and, for $d=1$, $\widehat{N}$ is the cardinality of $\widehat{\cal M}$ (cf. Remark 2.1f) ). Note that when $d\geq 2$, $\widehat{a}^{(d)}_{1-\alpha}$ has a typical structure appearing in confidence bands for probability density or regression functions. See, for example, the main result of Bickel and Rosenblatt (1973). When $d=1$, $\widehat{a}^{(d)}_{1-\alpha}$ corresponds to a quantile value of the maximum of a Gaussian mixture model because $\cal M$ is a collection of separated points under our assumptions. With this notation, our first confidence interval based on vertical variation, and using the de-biased estimator of the underlying density, is defined as

[TABLE]

Theorem 2.1

Suppose that (F1), (F2), (K), (A), (H1)0 and (H2)0 hold. Let $0<\alpha<1.$ If $\widehat{\mathscr{V}_{d-1}(\mathcal{M})}$ is a consistent estimator for $\mathscr{V}_{d-1}(\mathcal{M})$, then we have,

[TABLE]

An asymptotic confidence region for the superlevel set $\mathcal{L}$ is given by the following upper and lower bounds:

[TABLE]

and

[TABLE]

Corollary 2.1

Suppose that (F1), (F2), (K), (A), (H1)0 and (H2)0 hold. Also suppose $\beta_{n,h}^{(1),E}\rightarrow 0$ and $\beta_{n,l}^{(3),E}\rightarrow 0$ as $n\rightarrow\infty$. Let $0<\alpha<1.$ If $\widehat{\mathscr{V}_{d-1}(\mathcal{M})}$ converges to $\mathscr{V}_{d-1}(\mathcal{M})$ in probability, then we have

[TABLE]

Remark 2.2

a) Notice that the construction of $\widehat{C}_{n,1}$ involves the choice of two bandwidths, $h$ and $l$. This is the case for all the confidence regions considered in this paper that are based on $\widehat{f}^{\,bc}$.

b) Constructing $\widehat{C}_{n,1}(1-\alpha)$ with $\widehat{f}$ rather than $\widehat{f}^{\,bc}$ also results in an asymptotically valid confidence set for ${\cal M},$ provided undersmoothing is being used to handle the bias. In this case, a sufficient condition for consistency of the coverage probability is that the bandwidth satisfies $\beta_{n,h}^{(2),E}\rightarrow 0$. When this assumption holds, the stochastic term $\sup_{x\in\mathcal{M}}|\widehat{f}(x)-\mathbb{E}\widehat{f}(x)|$ asymptotically dominates the bias term $\sup_{x\in\mathcal{M}}|\mathbb{E}\widehat{f}(x)-f(x)|$, so that the latter can be ignored in the proof (see Hall, 1993).

c) We discuss two choices for the estimator $\widehat{\mathscr{V}_{d-1}(\mathcal{M})}$. Let $\lambda$ be the $d$-dimensional Lebesgue meaure, and let $P_{n}=n^{-1}\sum_{i=1}^{n}\delta_{X_{i}}$ denote the empirical probability measure, where $\delta_{x}$ denotes Dirac measure in $x$. Let $\mathscr{A}$ be the class of compact sets with a positive reach bounded away from zero such that $\mathcal{L}\in\mathscr{A}$. One option for the consistent estimator of $\mathscr{V}_{d-1}(\mathcal{M})$ is given by $\mathscr{V}_{d-1}(\partial\widehat{\mathcal{L}})$, where*

[TABLE]

The consistency of this estimator is shown in Proposition 3 in Cuevas et al. (2012). There, consistency is derived in terms of outer Minkowski content, which, under our assumptions, is equivalent to the consistency using Hausdorff measure (see Corollary 1 in Ambrosio et al., 2008). Efficiently computing $\widehat{\mathcal{L}}$ is challenging.

Another estimator for $\mathscr{V}_{d-1}(\mathcal{M})$ is given by $\mathscr{V}_{d-1}(\widehat{\mathcal{M}})$ with $\widehat{\mathcal{M}}=\widehat{f}^{-1}(c).$ The convergence rate and asymptotic normality of this estimator in the context of surface integral estimation are shown in Theorem 1 in Qiao (2018b), where additional assumptions are imposed, in particular on the speed of convergence of $h$.

d) The $(d-1)$-dimensional Hausdorff measure of $\mathcal{M}$ and its estimators come into play because (i) the distribution of the statistic related to the above confidence regions can be approximated by that of the extreme value of certain Gaussian random fields indexed by the level set; and (ii) the latter, in turn, is related to an integral over $\mathcal{M}$ with respect to the $(d-1)$-dimensional Hausdorff measure. In our case, the integrand of this integral turns out to be a constant - see Theorem 6.1, and thus we obtain the volume of the isosurface. In fact, it is well known that the probability of the extreme value of a locally stationary Gaussian random fields exceeding a large level is asymptotically proportional to the volume of the index set (locally speaking). See, for example, Chapter 2 of Piterbarg (1996). Surface integrals have appeared in the context of level-set estimation before in Cadre (2006). There, however, a first order asymptotics (consistency) is considered, using the set-theoretic measure of symmetric difference $d({\cal L},\widehat{\cal L})$. On a very heuristic level, the fact that a surface integral comes in here can be understood by the fact that the variance of the fluctuations of $d({\cal L},\widehat{\cal L})$ is of the order $a_{n}=1/nh^{d},$ which is inherited from the fluctuations of the density estimator, and by then approximating $d({\cal L},\widehat{\cal L})$ by a constant times {\rm vol}\big{(}f^{-1}[c+a_{n},c-a_{n}]\big{)}, the Lebesgue measure of $f^{-1}[c+a_{n},c-a_{n}].$ Lebesgue’s theorem gives that a^{-1}_{n}{\rm vol}\big{(}f^{-1}[c+a_{n},c-a_{n}]\big{)} converges to a surface integral over the boundary ${\cal M}$. In other words, the technical reason for this surface integral to appear in Cadre (2006) is different from why it appears in our context. **

2.2 Bootstrap confidence regions

Bootstrap confidence regions based on vertical variation of the kernel density estimate have been constructed in Mammen and Polonik (2013) and in Chen et al. (2017). They are based on a bootstrap approximation of quantiles of statistics of the form

[TABLE]

where $D_{n}$ is such that (asymptotically) ${\cal M}=f^{-1}(c)\subset D_{n}.$ The confidence regions considered in the literature differ in the choice of the set $D_{n}.$ While Mammen and Polonik (2013) propose to use $D_{n}=\{c-\epsilon_{n}\leq f(x)\leq c+\epsilon_{n}\}$ for some appropriate choice of $\epsilon_{n}$ tending to zero, as $n\to\infty$, Chen et al. simply use $D_{n}={\mathbb{R}}^{d}.$ Here we are using the smallest such set $D_{n}={\cal M}.$ Clearly, the statistics are stochastically ordered in terms of the size of the set $D_{n}$. Thus, the choice $D_{n}={\cal M}$ leads to the confidence set that is the smallest among the three. Of course the coverage still needs to be investigated. However, if the corresponding bootstrap approximations of the distributions of $T(D_{n})$ (and $T^{E}(D_{n})$) work similarly well, then using the statistic $T({\cal M})$ or $T^{E}({\cal M})$, respectively, can be expected to be a good choice for the construction the bootstrap based confidence sets based on the vertical variation.

Our construction is as follows. Let $X_{1}^{*},\ldots,X_{n}^{*}$ be a sample drawn from a kernel density estimator $\widehat{f}_{g}$ using $X_{1},\ldots,X_{n},$ and bandwidth $g>0$. Let $\widehat{f}^{*}(x)$ be the kernel density estimate using $X_{1}^{*},\ldots,X_{n}^{*}$ and bandwidth $h$. Let $\widehat{f}^{*,E}(x)=\mathbb{E}^{*}\widehat{f}^{*}(x)$, where we use $\mathbb{E}^{*}$ to indicate the expectation with respect to $\widehat{f}_{g}$. For $0<\alpha<1$, let $\widehat{c}_{1-\alpha}^{*,E}$ be the $(1-\alpha)$-quantile of the distribution of $\sup_{x\in\widehat{\mathcal{M}}}|\widehat{f}^{*}(x)-\widehat{f}^{*,E}(x)|$, and let $\widehat{c}_{1-\alpha}^{*}$ be the corresponding quantile of $\sup_{x\in\widehat{\mathcal{M}}}|\widehat{f}^{*}(x)-\widehat{f}_{g}(x)|$. We now define our bootstrap confidence regions for ${\cal M}$ and ${\cal M}^{E}$, respectively, as

[TABLE]

We also define the following sets to construct bootstrap confidence regions for ${\cal L}$ and ${\cal L}^{E}$, respectively.

[TABLE]

Below we show that this (and other) confidence region is asymptotically exact, and we derive rates of convergence for the coverage probability. These rates of convergence have a somewhat complex appearance, which we first explain from a high level perspective.

Structure of the rates of convergence for the coverage probabilities of bootstrap based confidence sets: The derivation of the following somewhat complex looking rates of convergence of the coverage probabilities are all based on Lemma 8.1, which is a slight reformulation of a result of Mammen and Polonik (2013). Based on this result, the rates are of the form O\big{(}\sqrt{\delta_{n}}+\tau_{n}), where $\delta_{n}$ and $\tau_{n}$ are derived as follows: Let Z_{n}=\sup_{x\in A}\big{|}\widehat{g}(x)-\overline{g}(x)\big{|}, and let $Z_{n}^{*}$ be a bootstrap version, where $\widehat{g}(x)$ is one of the density estimators considered, and $\overline{g}(x)$ is some centering; the set $A$ in the supremum is either $\cal M$ or ${\cal M}^{E}.$ Then, we derive $\delta_{n}$ and $\tau_{n}$ by showing that for some sequence of positive real numbers $\gamma_{n}$, we have P\big{(}\big{|}Z_{n}-Z_{n}^{*}\big{|}>\gamma_{n}\big{)}\leq\delta_{n} and \sup_{t\in{\mathbb{R}}}P\big{(}Z_{n}\in[t,t+\gamma_{n})\big{)}\leq\tau_{n}. Both of $\delta_{n}$ and $\tau_{n}$ are themselves comprised of a sum of various terms. In fact, in our applications of this result, $\gamma_{n}$ is chosen such that $\sqrt{\delta_{n}}$ is negligible, and we have $\tau_{n}=\Psi_{n}(\gamma_{n})$, with

[TABLE]

That, in our case, $\tau_{n}$ has this particular form follows from a result by Neumann (1998). The rates $\gamma_{n}$ that make $\sqrt{\delta_{n}}$ negligible will follow, respectively, from strong approximation results in Neumann (1998) and some modification of Neumann’s result given in Mammen and Polonik (2013).

Note further that in each of the construction approaches considered below, we obtain the same rates of convergence of the coverage probabilities for both confidence regions for ${\cal M}$ and confidence regions for ${\cal L}$ (and similarly for ${\cal M}^{E}$ and ${\cal L}^{E}$). The reason for this is as follows. Let $p_{n,\cal M}$ and $p_{n,\cal L}$ denote these coverage probabilities based on one approach (e.g. the left-hand sides of (2.11) and (2.12), respectively. By construction, we have $p_{n,\cal L}\leq p_{n,\cal M}$. More precisely, it is shown in the proof of Corollary 2.1 that $p_{n,\cal L}=p_{n,\cal M}\cdot q_{n}$, where $0\leq 1-q_{n}=O(\delta_{n})$ with $\delta_{n}=o(1).$ Now, if $p_{n,\cal M}=1-\alpha+O(\tau_{n}),$ then we obtain p_{n,\cal L}=(1-\alpha+O(\tau_{n}))(1-O(\delta_{n}))=(1-\alpha)+O\big{(}\max(\tau_{n},\delta_{n})\big{)}, and under our respective assumptions, O\big{(}\max(\tau_{n},\delta_{n})\big{)}=O(\tau_{n}), showing that both $p_{n,\cal M}$ and $p_{n,\cal L}$ converge to $(1-\alpha)$ at the same speed.

Theorem 2.2

Suppose that (F1), (F2) and (K) hold. Let $0<\alpha<1,$ and

[TABLE]

(a) If $\Psi_{n}(\gamma_{n}^{E})=o(1),$ then we have

[TABLE]

(b) If $\Psi_{n}(\gamma_{n})=o(1)$ as $n\rightarrow\infty,$ then we have

[TABLE]

Remark 2.3

a) The set $\widehat{C}_{n,2}^{*,E}$ can be also used as a confidence region for $\mathcal{M}$ (not just for $\mathcal{M}^{E}$) if the bandwidth is chosen to be of smaller order than the optimal bandwidth, to make the bias negligible (undersmoothing). In practice, the choice of undersmoothing bandwidth might be difficult to determine. We do not pursue the theoretical justification for $\widehat{C}_{n,2}^{*,E}$ as a confidence region for $\mathcal{M}$, the numerical performance of which, however, is shown in the simulation section.

b) The quantity $\beta_{n,g}^{(2)}$ appears in $\gamma_{n}$, because the second partial derivatives appear in the bias of density estimation and need to be estimated in our proof. Note that we are not using the de-biased density estimator in this theorem.

c) For $d\geq 2$, if we choose both $g$ and $h$ to be of the standard optimal rates, i.e. $g=h={\rm const.}\,n^{-1/(d+4)}$, then

[TABLE]

The rate in Chen et al. (2017) is given by $(nh^{d})^{-1/8}(\log{n})^{7/8}$. When $h$ is chosen as the standard optimal bandwidth for density estimation, i.e. $h$ is of the exact order $n^{-1/(d+4)}$, the rate given in Chen et al. (2017) becomes $n^{-1/(2d+8)}(\log{n})^{7/8}$, which is slower than $\rho_{n}^{E}$. Note that even if the rate in Chen et al. (2017) is $(nh^{d})^{-1/6}(\log{n})^{7/6}$ as claimed in this paper, it is still slower than $\rho_{n}^{E}$ when using the optimal bandwidth.

For $d\geq 2$, if we again choose $h=O(n^{-1/(d+4)})$, then

[TABLE]

With $g=O(n^{-1/(d+6)})$, $\rho_{n}$ is of the order of $n^{-1/(d+6)}\log{n}$. **

2.2.1 A bootstrap confidence region based on explicit bias correction

Introducing a new bandwidth $g$ and bootstrapping from $\hat{f}_{g},$ as we did above, can be viewed a method of bias correction (see page 208, Hall 1992). This allows us to construct a confidence region for ${\cal M}$ (rather than just for ${\cal M}^{E}$). We can also construct a confidence region for $\mathcal{M}$ using an explicit bias correction by

[TABLE]

For confidence regions for ${\cal L}$, we define

[TABLE]

We have the following result:

Theorem 2.3

Suppose that (F1), (F2), and (K) hold. Let $0<\alpha<1.$ Let

[TABLE]

If $\beta_{n,l}^{(2)}=o(1)$ and \Psi_{n}\big{(}\gamma_{n}^{bc}\big{)}=o(1) as $n\rightarrow\infty$, then we have

[TABLE]

If we further assume $\beta_{n,l}^{(3)}=o(1)$ as $n\rightarrow\infty$, then

[TABLE]

Remark 2.4

The constructions of both $\widehat{C}_{n,2}^{*,E}$ and $\widehat{C}_{n,3}^{*}$ are using a quantile, $\widehat{c}^{*,E}_{1-\alpha},$ that is ignoring the bias. Nevertheless, $\widehat{C}_{n,2}^{*,E}$ gives a confidence region for the smoothed isosurface ${\cal M}^{E}$, while $\widehat{C}_{n,3}^{*}$ is a confidence region for ${\cal M}$. This is so, because one of them, $\widehat{C}_{n,3}^{*},$ is based on the de-biased estimator, while $\widehat{C}_{n,2}^{*,E}$ is not. Heuristically, this can be understood by writing $\widehat{f}^{\,bc}(x)-f(x)=\widehat{f}(x)-\mathbb{E}\widehat{f}(x)+\beta(x)-\widehat{\beta}(x).$ One can see that the bias correction in the density will adjust for the bias that is present in the quantile $\widehat{c}^{*,E}_{1-\alpha}.$

For $d\geq 2$, if we choose the optimal bandwidth $g=h=O(n^{-1/(d+4)})$, and $l=O(n^{-1/(d+8)})$ (which is the order of the optimal bandwidth for estimating the second derivatives), then

[TABLE]

Compared to the rate $\rho_{n}$ given in Remark 2.3 this rate is faster. However, the fact that the construction of $\widehat{C}_{n,3}^{*}$ involves the choice of three bandwidths might be a challenge in practice. **

3 Confidence regions based on horizontal variation

Various confidence regions for ${\cal M},{\cal M}^{E},{\cal L}$ and ${\cal L}^{E}$ based on horizontal variation will be derived in this section. The geometric link between horizontal and vertical variation based confidence regions is, for obvious reasons, given by the gradient. Let $x\in{\cal M}$ and $\widehat{x}\in\widehat{\cal M}$, such that $\widehat{x}$ is close to $x$, then we obviously have

[TABLE]

and when $\widehat{x}$ is chosen such that $x-\widehat{x}$ is approximately perpendicular to $\widehat{\cal M}$,

[TABLE]

Different ways of choosing $\widehat{x}$ will give rise to different types of horizontal variation based confidence regions. For example, $\widehat{x}$ can be a projection point of $x$ onto $\widehat{\cal M}$ or it can be chosen such that there exists a gradient integral curve connecting $x$ and $\widehat{x}$.

While the above confidence regions based on vertical variation are using approximations to quantiles of the distribution of $\sup_{x\in{\cal M}}|\widehat{f}(x)-f(x)|$, the horizontal variation based confidence regions will use estimated quantiles of the quantity \sup_{x\in{\cal M}}\big{|}\nabla\widehat{f}(\widehat{x})^{T}\big{(}x-\widehat{x}\big{)}\big{|}. The latter methods are not purely horizontal variation based, as they involve the adjustment of $\widehat{x}-x$ by the gradient of $f$. Nevertheless, since they are explicitly using the differences $x-\widehat{x}$, we still call them methods based on horizontal variation.

Confidence regions based on horizontal variation (without estimating the gradient), have been constructed in Chen et al. (2017) based on the Hausdorff distance. The approaches considered in our work are asymptotically equivalent to the method proposed in Chen et al. (2017), but in practice the confidence sets are different. The different constructions also provide additional insight into the underlying geometry.

In the following we introduce novel horizontal variation based methods, based on estimating quantiles via both the asymptotic distribution and the bootstrap. Instead of using standard bootstrap as in Chen et al. (2017), we adopt smoothed bootstrap. We would like to note that rather than simply providing alternative methods for confidence regions, the constructions and the discussion of the performance of the resulting confidence regions are meant to provide insight into the effects of geometric aspects of the underlying probability density on the performance of confidence regions.

3.1 Methods based on integral curves

The approach discussed here is based on the construction of $\widehat{x}$ (cf. (3.1)) using integral curves driven by the (scaled) gradient field and their relation to level set. This relation will be discussed first.

Integral curves and level sets: For any $x_{0}\in\mathbb{R}^{d}$, let $\{\mathfrak{X}_{x_{0}}(t),t\in\mathbb{R}\}$ be the integral curve driven by the scaled gradient of $f,$ starting from $x_{0}$, defined by the equation,

[TABLE]

where we assume that $\|\nabla f(\mathfrak{X}_{x_{0}}(t))\|\neq 0$ (cf. Assumption (F2)). For $d=1$, $\nabla f$ is understood to mean $f^{\prime}$. The reason for choosing the scaled gradient $\frac{\nabla f}{\|\nabla f\|^{2}}$ as a driving vector field, rather than $\nabla f$ itself, is based on the following convenient property. Suppose that $t>0$, and set $I_{0}(t)=[0,t]$. Then, by the fundamental theorem of calculus for line integrals, we have for any $x_{0}\in\cal M$ and $t\in\mathbb{R}$ such that $\|\nabla f(x)\|\neq 0$ for all $x\in\{s:c\leq f(s)\leq c+t\},$

[TABLE]

where the right hand side of (3.3) is a line integral over the trajectory $\{\mathfrak{X}_{x_{0}}(s):s\in I_{0}(t)\}$, and “$\cdot$” represents dot product between vectors. By adopting the current scaled vector field, we can make sure the height increase (or decrease) in the density is exactly the amount of “time” needed to travel. In other words, if two particles start from any two points $x_{1},x_{2}\in\mathcal{M}$, then, after following the integral curves $\mathfrak{X}_{x_{1}}(\cdot)$ and $\mathfrak{X}_{x_{2}}(\cdot)$ respectively for time $t$ (which can be negative), both of these two particles will arrive at $\mathcal{M}_{c+t}$. The same holds for $t\leq 0$, by using the convention to start integration at [math].

Let the plug-in estimators of $\nabla f(x)$ and $\mathfrak{X}_{x_{0}}(t)$ based on the kernel density estimate be denoted by $\nabla\widehat{f}(x)$ and $\widehat{\mathfrak{X}}_{x_{0}}(t)$, respectively, where the latter is the solution of the differential equation

[TABLE]

We use the notation $\widehat{\mathfrak{X}}^{bc}_{x_{0}}(t)$ to denote the integral curve as in (3.5), but with $\widehat{f}$ replaced by the bias-corrected version $\widehat{f}^{\,bc}$.

3.1.1 Confidence regions for ${\cal M}$ and ${\cal L}$ using local adjustment by the gradient

Using the de-biased estimator $\widehat{f}^{\,bc}$ along with the integral curve $\widehat{\mathfrak{X}}^{bc}_{x}(t)$, we now present the construction of two confidence regions, one based on asymptotic distribution theory, and the other based on the bootstrap. The latter has a faster rate of approximation of the coverage probability.

For $x\in{\mathbb{R}}^{d},$ define $\widehat{\theta}^{\,bc}_{x}\in{\mathbb{R}}$ through the property

[TABLE]

For large sample size the existence and uniqueness of $\widehat{\theta}^{\,bc}_{x}$ for $x\in\mathcal{M}$ are proved in Lemma 7.1. For finite sample, in case the solution to (3.6) is not unique, we take $\widehat{\theta}^{\,bc}_{x}$ as the infimum of the set of solutions; and whenever there is no solution to (3.6), we set $\widehat{\theta}^{\,bc}_{x}$ to be the smallest value of $\mathop{\rm argmin}_{\theta}|\widehat{f}^{\,bc}(\widehat{\mathfrak{X}}^{bc}_{x}(\theta))-c|$. Noticing that $\widehat{\mathfrak{X}}^{bc}_{x}(\widehat{\theta}^{\,bc}_{x})\in\widehat{\cal M}^{\,bc}=(\widehat{f}^{\,bc})^{-1}(c)$, we now set $\widehat{x}$ in (3.1) as $\widehat{x}=\widehat{\mathfrak{X}}^{bc}_{x}(\widehat{\theta}^{\,bc}_{x})$.

Letting $0<\alpha<1,$ and recalling the definition of $\widehat{a}^{(d)}_{1-\alpha}$ given in (2.3), we define

[TABLE]

It can be shown (see Lemma 7.1) that $x\mapsto\widehat{\mathfrak{X}}^{bc}_{x}(\widehat{\theta}^{\,bc}_{x})$ defines a bijective map between $\mathcal{M}$ and $\widehat{\mathcal{M}}^{bc}$ when the sample size is large enough. Thus, for large sample size, we can equivalently write

[TABLE]

This also indicates an algorithm: For a dense enough subset of values $z\in\widehat{\cal M}^{\,bc}$, run the integral curve $\widehat{\mathfrak{X}}^{\,bc}_{z}(t)$, and check whether the condition in the definition of the confidence region holds. A bootstrap version of this confidence region is given by

[TABLE]

where $\widehat{c}^{*,E}_{1-\alpha}$ is as in (2.2).

Let $\widehat{\mathcal{L}}^{bc}=\{x\in\mathbb{R}^{d}:\;\widehat{f}^{bc}(x)\geq c\}$ and define

[TABLE]

Theorem 3.1

Part 1.* Let $0<\alpha<1.$ Suppose that (F1), (F2), (K), (A), (H1)2 and (H2)2 hold. Then we have, as $n\to\infty$,*

[TABLE]

Part 2.* Let $0<\alpha<1,$ and suppose that (F1), (F2), (K), (H1)2 and (H2)2 hold. Let $\gamma_{n}^{bc}$ as in (2.13). If $\Psi_{n}(\gamma_{n}^{bc})=o(1)$ as $n\rightarrow\infty$, then we have*

[TABLE]

Notice that we do not have rate of convergence for the coverage probability of the confidence regions based on the large sample theory. At this point it is an open question whether, and if yes, how, to derive such rates of approximations. Since we use the approximation of the extreme value distribution of Gaussian fields, it is expected to be a slow rate $(\log{n})^{-1}$, following a similar argument given in Hall (1979).

3.1.2 Confidence regions for ${\cal M}^{E}$ and ${\cal L}^{E}$ adjusted by gradient

Similar to the above, confidence regions for ${\cal{M}}^{E}$ and ${\cal{L}}^{E}$ can be constructed and analysed. For instance, a bootstrap confidence region for ${\cal{M}}^{E}$ is given by

[TABLE]

The lower and upper “bounds” of the confidence region for ${\cal L}^{E}$ can be constructed by taking set difference and union between $\widehat{\cal L}$ and $\widehat{C}^{*,E}_{n,4}(1-\alpha)$, respectively. Note that the construction uses estimators based on $\widehat{f}$, instead of $\widehat{f}^{\,bc}$, similar to $\widehat{C}^{*,E}_{n,2}(1-\alpha)$ in (2.2). It can be shown that their rates of convergence for the coverage probability are O\big{(}\Psi(\gamma_{n}^{E})\big{)} if we assume that (F1), (F2), (K), (H1)2 and (H2)2 hold and $\Psi(\gamma_{n}^{E}))=o(1)$. The proof follows the same arguments given in the proof of Theorem 3.1. Details are omitted.

3.1.3 Confidence regions not locally adjusted by the gradient

The confidence regions constructed here are related to the confidence regions constructed in the previous subsection, but in contrast to them, here the dependence on the estimated gradient is more indirect through the construction of the integral curve. As a result, the width of the confidence regions only depends on the length of the integral curve, and it is not locally adjusted by the gradient. The construction is as follows. Recall that $\widehat{f}^{*,E}(x)={\mathbb{E}}^{*}\widehat{f}^{*}(x),$ and let

[TABLE]

Let $\widehat{\mathfrak{X}}_{x}^{*}$ be the integral curve driven by $\nabla\widehat{f}^{*}$, and let $\widehat{\theta}_{x}^{*}$ be the first time $t$ at which $\widehat{\mathfrak{X}}_{x}^{*}(t)$ hits $\mathcal{M}^{*}$. In the case of large sample size the existence and uniqueness of $\widehat{\theta}_{x}^{*}$ follow the same argument as for $\widehat{\theta}_{x}$. For $0<\alpha<1,$ let $\widehat{d}_{1-\alpha}^{*,E}$ be the quantile of order $(1-\alpha)$ for $\sup_{x\in\widehat{\mathcal{M}}^{*,E}}|\widehat{\mathfrak{X}}_{x}^{*}(\widehat{\theta}_{x}^{*})-x|,$ and define the set

[TABLE]

Also define

[TABLE]

Let $\zeta_{n}=\beta_{n,g}^{(1)}\log{n}+\beta_{n,g}^{(0)}\left[h^{-1}\log{n}+\sqrt{nh^{d+4}\log{n}}\right]$. Now we have the following result:

Theorem 3.2

Suppose that (F1), (F2), (K), (H1)2 and (H2)2 hold. Then, if $\Psi_{n}(\gamma_{n}^{bc})+\zeta_{n}=o(1)$ as $n\rightarrow\infty,$ where $\gamma_{n}^{bc}$ is as in (2.13), we have as $n\rightarrow\infty$,

[TABLE]

Remark 3.1

Note that if we further assume that $g$ and $h$ are of the same rate, i.e., there exist $0<C_{1},C_{2}<\infty$ such that $C_{1}<h/g<C_{2}$ as $n\rightarrow\infty$, then $\zeta_{n}$ can be absorbed into $\Psi_{n}(\gamma_{n}^{bc})$ in the above results.

Similar to subsection 3.1.2, a confidence region for ${\cal M}^{E}$ can be constructed by

[TABLE]

and the confidence region for ${\cal L}^{E}$ has ‘upper and lower boundaries’ given by

[TABLE]

3.2 Horizontal variation based methods not based on integral curves

It follows from the Tubular Neighborhood Theorem (e.g. see Theorem 11.4 of Bredon, 1993) that for all $x\in\mathcal{M}$ and $n$ large enough, there exist unique $X_{x}\in\widehat{\mathcal{M}}^{bc}$ and $s_{x}\in\mathbb{R}$ such that $x=X_{x}+s_{x}\nabla\widehat{f}^{bc}(X_{x})$. Similarly, for all $x\in\widehat{\mathcal{M}}^{*,E}$ and large sample size, we can find $X_{x}^{*}\in\widehat{\mathcal{M}}^{*}$ and $s_{x}^{*}$ be such that $x=X_{x}^{*}+s_{x}^{*}\nabla\widehat{f}^{*}(X_{x}^{*})$.

Now let $\widehat{b}_{1-\alpha}^{*,E}$ be the quantile of $(1-\alpha)$ for $\sup_{x\in\widehat{\mathcal{M}}^{*,E}}\|X_{x}^{*}-x\|$, and for $0<\alpha<1$ define

[TABLE]

Also define

[TABLE]

Under suitable regularity conditions $\widehat{C}^{*}_{n,6}(1-\alpha)$ is a confidence region for ${\cal M}$ of asymptotic coverage level $1-\alpha$, and $\widehat{C}_{n,6}^{*,-}(1-\alpha)$ and $\widehat{C}_{n,6}^{*,+}(1-\alpha)$ give the upper and lower “bounds” of an asymptotic $(1-\alpha)$ confidence region for $\mathcal{L}$. The convergence rate of the coverage probability of $\widehat{C}^{*}_{n,6}(1-\alpha)$ for $\mathcal{M}$ as well as $\{\widehat{C}_{n,6}^{*,+}(1-\alpha),\widehat{C}_{n,6}^{*,-}(1-\alpha)\}$ for $\mathcal{L}$ can be derived in the way similar to the proof of Theorem 3.2. However, we do not further pursue it here. While the geometric construction of the confidence region $\widehat{C}^{*}_{n,6}(1-\alpha)$ is essentially the same as the one constructed in Chen et al. (2017), our derivation via the tubular neighborhood theorem provides a slightly different angle to the construction. A similar confidence region for ${\cal M}^{E}$ is given by

[TABLE]

For ${\cal L}^{E}$, we can use the following upper and lower “bounds”:

[TABLE]

Rates of convergence of the coverage probabilities for these confidence regions can be derived by using similar ideas as above. No further details are given.

4 Performance of confidence regions and geometry

The above discusses the large sample behavior of various confidence regions for level sets. Bootstrap based methods show a faster rate of convergence of their coverage probability to the nominal level than the methods based on asymptotic distribution theory, which is not a surprise. A more detailed comparison based on the theoretical developments is not entirely straightforward, because different confidence sets depend on a different number of bandwidths to be chosen, requiring different assumptions, etc. However, for finite samples, certain relations between the geometry of the underlying density and the performance of the different types of confidence regions give some interesting insight. This will be discussed now.

1. By construction, most of the confidence regions for ${\cal M}$ or ${\cal M}^{E}$ based on horizontal variation constitute a band (or a tube) of constant width about the estimated target isosurface. The width of the tube depends on the global behavior of the density in a neighborhood about the targeted isosurface. Since horizontal variation based methods are essentially based on the worst case behavior (similar to the supremum distance), it can be expected that for densities for which the norm of the gradient varies a lot along the isofurface, this confidence band tends to be unnecessarily wide. This is also illustrated in our simulation study in Section 5.

In contrast to that behavior of the horizontal confidence regions, the width of the vertical variation based confidence regions has local adaptivity. Essentially, their width at a given point of the isosurface is inversely proportional to the norm of the gradient, and thus they are containing additional information about the geometry of the density.

2. Our theoretical assumptions restrict the level $c$ to be strictly larger than zero, which means that in local neighborhoods about the corresponding isosurface the density is bounded away from zero. However, a relatively small level $c$ can still provide problems in finite samples. For instance, vertical variation based methods of the form $\widehat{f}[c-\widehat{a}_{n},c+\widehat{a}_{n}]$ (or the ones using the bias-corrected density estimator) might be very large in volume, as the lower bound $c-\widehat{a}_{n}$ might be less than zero, meaning that the outer confidence region only ends at the support of the density estimator used. Nevertheless, the probability content carried by the confidence regions might still be small. (Asymptotically, this problem of course disappears simply because $\widehat{a}_{n}$ converges to zero.) In such situations, the volume of horizontal variation based methods tend to be of smaller volume than the ones of the vertical based methods, but the probability mass carried by them might nevertheless be larger. This can be seen in the simulation results presented in Table 1 when inspecting the column corresponding to Case 3.

3. Another interesting scenario corresponds to levels $c$ that are close to critical values of the density. Similar to the previous item, this problem does not appear in a large sample scenario, because our assumptions require the gradient along the iso-surface to be bounded away from zero. For finite samples, however, we observe the following interesting geometric challenge.

Suppose that $c$ is only slightly smaller than a level at which the true density has a local maximum that is not a global maximum, and let $x_{0}$ denote the point at which this local maximum is attained. Then, under our regularity assumptions, there will be a neighborhood of $x_{0}$ that is part of the superlevel set ${\cal L}$. However, the value of the density estimates $\widehat{f}(x_{0})$ or $\widehat{f}^{\,bc}(x_{0})$ might, due to random fluctuation, not exceed the value $c$, so that the estimated superlevel set does not contain a neighborhood about $x_{0}$. Since the confidence regions based on horizontal variation are built around the contours of these estimated superlevel sets, they might miss such areas - and note that these areas then could in fact lie far from the confidence region. In such cases the missed target regions tend to be small in size. Nevertheless, the topology (homology) of the confidence regions will in general be different from the one of the target contour, and the normal compatibility assumption (e.g. see Chazal et al. 2007) used in Chen et al. (2017) for the construction of horizontal variation based confidence regions will be violated in such cases.

Horizontal variation based confidence regions will not perform well in this scenario. Observe that the quantiles used in their construction are essentially based on the maximum distance of the estimated and the true contours. This distance tending to be large means that these quantiles will tend to become large, leading to wide confidence tubes. By contrast, the vertical distance is less impacted by the different topology of the estimated and true contours, and hence the vertical variation based confidence regions suffer less from having very large volume in this scenario.

A scenario as the one just discussed is included in the simulation study presented below. See Table 1, Case 4, where one can see that the Hausdorff based methods tend to be quite large. A similar remark applies to the integral curve based methods, such as $\widehat{C}_{n,4}(1-\alpha)$, where the indicated problem is expressed by the non-existence of $\widehat{\theta}_{x}^{\,bc}$ for a non-negligible set of starting values $x$. A similar discussion applies when the level $c$ is narrowly above a critical level, illustrated in Figure 1.

4. Confidence regions based on estimated integral curves might, for finite sample size, suffer from the local geometry of the kernel density estimator, as is illustrated in Figure 2. The two panels in this figure show a scenario (for finite sample size) that violates the bijective condition for horizontal methods, which is shown to hold for large samples (see Lemma 7.1).

The violation is due to small local minimum of the kernel density estimator, and this local geometric property then ‘diverts’ the integral curves from their expected path. This also gives rise to numerical challenges. In the situation shown in the figure, a sample of size 200 was drawn from a density function in (5.1) below with $a=2,$ and then a bootstrap sample was drawn. The focus was on the level set corresponding to $p=0.95$. In the left panel, the pink curves are contour lines of $\widehat{f}$; the red curve is $\widehat{\mathcal{M}}$; the black curve is $\mathcal{M}$; the blue curves are trajectories of integral curves driven by the gradients of $\widehat{f}$. In the right panel, the pink curves are contour lines of $\widehat{f}^{*}$; the red curve is still $\widehat{\mathcal{M}}$; the green curve is $\mathcal{M}^{*}$; the blue curves are trajectories of integral curves driven by the gradients of $\widehat{f}^{*}$. Notice that the trajectories fail to define a bijective mapping between $\widehat{\mathcal{M}}$ and $\widehat{\mathcal{M}}^{*}$ around (-1.5, -1.5) due to the existence of a local minimum of $\widehat{f}^{*}$.

5 Simulations

This simulation study compares six (one large sample based, and five bootstrap based) confidence regions for ${\cal M}$ and ${\cal M}^{E}$ in terms of coverage probability and volume, thereby

• •

considering different levels $c$ (low, high, close to critical levels), and

• •

comparing vertical variation based and horizontal variation based methods.

All the horizontal variation based methods are more computationally involved than the vertical methods. A preliminary simulation study showed our horizontal variation based methods to behave similarly to the Hausdorff-distance based approach of Chen et al. (2017). Therefore we here only use the latter to represent the horizontal variation based methods. With $d(x,\widehat{\cal M})=\inf_{s\in\widehat{\cal M}}\|x-s\|,$ these confidence regions have the form \big{\{}x\in{\mathbb{R}}^{d}:d(x,\widehat{\cal M})\leq\widehat{e}^{*}_{1-\alpha}\big{\}}, where $\widehat{e}^{*}_{1-\alpha}$ is a bootstrap based estimate of the $(1-\alpha)$-quantile of $d(x,\widehat{\cal M}).$ Recall that the Hausdorff distance between ${\cal M}$ and $\widehat{\cal M}$ is given by

[TABLE]

and thus, $\sup_{x\in{\cal M}}d(x,\widehat{\cal M})$ is ‘one part’ of the Hausdorff distance. However, if ${\cal M}$ and ${\widehat{\cal M}}$ are normal compatible, then $\sup_{x\in{\cal M}}d(x,\widehat{\cal M})=\sup_{s\in\widehat{\cal M}}d(s,{\cal M})$, and

[TABLE]

Chen et al. (2017) are using this approach with ${\cal M}$ replaced by ${\cal M}^{E}$, and they show that, under certain regularity assumptions, normal compatibility of ${\cal M}^{E}$ and ${\widehat{\cal M}}$ holds asymptotically with probability tending to one.

One of the models used in our simulations is the bivariate normal with

[TABLE]

where the contours of the density function are ellipses with $a$ defining their eccentricity. Then, for $c=f(x_{0},y_{0};a)$ with $a^{2}x_{0}^{2}+y_{0}^{2}/a^{2}=r_{0}^{2}$ for some $0<r_{0}<\infty$, the probability over the superlevel set $\{(x,y):\;\;f(x,y;a)\geq c\}$ is

[TABLE]

We choose $p=50\%$ and $p=95\%$. Our second model is a mixture of normal distributions of the form

[TABLE]

This density has two modes with corresponding heights 0.065 and 0.11, respectively. In our study we chose the level $c=0.048$, which lies slightly below the lower local maximum of the mixture of normals (cf. 3. in Section 4). We consider the following 4 cases:

Case 1: density in (5.1) with $a=1$ and $p=0.5$,

Case 2: density in (5.1) with $a=2$ and $p=0.5$,

Case 3: density in (5.1) with $a=1$ and $p=0.95$,

Case 4: density in (5.2) with $c=0.048$.

As a kernel we choose the form

[TABLE]

We ran the simulation for 400 times. In each iteration, a sample of size $n$ was randomly drawn from the given distribution and then a bootstrap procedure based on 250 bootstrap re-samplings was performed to create the confidence regions using the following methods:

(H) Hausdorff-distance based approach of Chen et al. (2017) for the smoothed level set;

(V.e)  vertical variation based confidence region $\widehat{C}_{n,2}^{*,E}(1-\alpha)$ for the smoothed level set;

(V)  vertical variation based confidence region $\widehat{C}_{n,2}^{*}(1-\alpha)$ for the true level set;

(V.bc)  vertical variation based confidence region with bias correction $\widehat{C}_{n,3}^{*}(1-\alpha)$ for the true level set;

(V.us)  vertical variation based confidence region with undersmoothing $\widehat{C}_{n,2}^{*,E}(1-\alpha)$ for the true level set (see Remark 2.3);

(V.ls)  vertical variation based large sample confidence region $\widehat{C}_{n,1}(1-\alpha)$ for the true level set.

The bandwidths involved in the construction of these confidence regions are selected using the direct plug-in method. In particular, we use the plug-in optimal bandwidth for kernel density estimation as $h$ and $g$, while using the plug-in optimal bandwidth for the second derivative estimation as $l$. In fact, we choose different bandwidths for each of the two dimensions. For (V.us), we used 70% of the optimal plug-in bandwidths.

The confidence level was set to be 90%. With the 400 runs, we calculated the coverage probabilities (C.P.) of these confidence regions as well as their average Lebesgue measures ($\bar{\lambda}$) and average probability measures ($\bar{P}$). Note that in the general form $\widehat{f}[c-\widehat{a},c+\widehat{a}]$ or $\widehat{f}^{bc}[c-\widehat{a},c+\widehat{a}]$ of the confidence regions based on vertical variation, sometimes $c-\widehat{a}<0$ for Case 3. Since the kernel function $K$ we used has bounded support, so do $\widehat{f}$ and $\widehat{f}^{bc}$. The outer boundary of the confidence regions based on vertical variation for the level sets is in fact the support of the density estimator (cf. 2. of Section 4). For numerical reasons, we used $\widehat{f}[\max(c-\widehat{a},\omega),c+\widehat{a}]$ or $\widehat{f}^{bc}[\max(c-\widehat{a},\omega),c+\widehat{a}]$ as the confidence regions, where we took $\omega=10^{-6}$.

In Table 1, (H) has to be compared with (V.e) since these methods are both targeting the smoothed level sets $\mathcal{M}^{E}$. Overall it is clear that the vertical method (V.e) outperforms the horizontal method (H). Detailed discussions have been given in Section 4.

We only include results for the large sample confidence regions (V.ls) with $n=5000$ and $n=50000$. This is because the formula in Theorem 2.1 requires $h<1$, which cannot be satisfied when $n$ is small in our examples. Overall the large sample confidence regions have conservative coverage probabilities in our examples. This is not surprising because it is well-known that the convergence rate of the coverage probability is slow for such large sample confidence regions. However, their volumes are comparable to those of the bootstrap confidence regions when the sample size is large, which makes (V.ls) a competitive option considering its computation does not require bootstrap.

Comparing (V), (V.bc) and (V.us), which are all confidence regions for the true level sets $\mathcal{M}$, it is apparent that (V) performs best. Convergence of the coverage probability of bias correction method (V.bc) is the slowest. This seems to be caused by slow convergence in the estimation of the second order derivatives in the bias correction. The confidence regions based on the undersmoothing method (V.us) have the largest volume. This is because the variance becomes large when the selected bandwidth is small.

6 Conclusion

We have constructed and analyzed various confidence regions for density superlevel sets and density isosurfaces based on plug-in estimates using kernel density estimation. The analysis is done in terms of large sample theory and also in the finite sample setting using simulations. The geometry underlying the construction of the different types of confidence regions is discussed. Geometric considerations also play a role in the interpretation of the finite sample behavior of the confidence regions. Overall, vertical variation based confidence intervals appear to have an edge over the horizontal methods.

The kernel estimator used in our investigations can of course be replaced by other (non-parametric) density estimators. For such modifications, the large sample behavior of the corresponding coverage probabilities might be analyzed using a similar Ansatz as in this work (see discussion of “Structure of the rates of convergence for the coverage probabilities of bootstrap based confidence sets” given in Section 2.2). This requires the investigation of all the relevant properties needed for our approach to go through.

There are various open questions related to the construction of confidence regions for density level sets. For instance, what can be said about optimality of the rates of convergence of the coverage probabilities? (Thanks to the referee for asking this question). Besides some classical work by Hall and Jing (1995), the only other work related to this question we are aware of is Calonico et al. (2018b). The role of bias correction in this context might be explored as well. For some recent work on bias correction see Chen (2017), and Calonico et al. (2018a). While we have been concentrating on density level sets, a similar approach might work for level sets of other functions, such as regression level sets, for instance. Level sets also play an integral role in the context of topological data analysis (persistent homology). Similar to our vertical variation based upper and lower confidence sets, Bobrowski et al. (2017) are using such upper and lower approximations of the level set to construct estimates for the topology of a single density level set. It might be worthwhile to explore this connection in more detail.

7 Proofs

7.1 Proof of Theorem 2.1

A key ingredient to the proof of Theorem 2.1, is the following special case of the main theorem in Qiao and Polonik (2018):

Theorem 7.1

Let $\mathcal{H}\subset\mathbb{R}^{d}$ ($d\geq 2$) be a compact set. Let $Z_{h}(x),x\in\mathcal{H}$, $0<h\leq 1$ be a sequence of centered Gaussian fields with covariance

[TABLE]

uniformly in $h\in(0,1]$ and $x\in\mathcal{H}$ as $\Delta x/h\rightarrow 0$, where $D$ is $d\times d$ positive definite matrix. Let $r<d$ and $\mathcal{M}\subset\mathcal{H}$ be a $r$-dimensional compact Riemannian manifold with reach $\Delta(\mathcal{M})>0$. For any $\delta>0$, define

[TABLE]

Suppose for any $\delta>0$, there exists a positive number $\eta$ such that

[TABLE]

In addition, assume that there exists $\eta>0$ and $\delta_{0}$, such that, for any $\delta>\delta_{0}$, we have

[TABLE]

For any fixed $z$, define

[TABLE]

where $M_{s}$ is a $d\times r$ matrix with orthonormal columns spanning ${\cal T}_{s}\mathcal{M}.$ Then

[TABLE]

This result will play a key role in the proof of Theorem 2.1, which is presented now. First we are going to prove (2.4) for $d\geq 2$. Recall $\|K\|_{2}^{2}=\int K^{2}(u)du$. Let

[TABLE]

Let

[TABLE]

To prove (2.4), using Slutsky’s Theorem, it suffices to show

[TABLE]

because we have

[TABLE]

To prove (7.5) we show the following two properties:

[TABLE]

and

[TABLE]

Using the uniform convergence rates for kernel density derivatives (see Lemma 3 in Arias-Castro et al. 2016 ) we obtain

[TABLE]

Property (7.8) now follows by using assumptions (H1)0 and (H2)0. Next we will show (7.7). With \mathcal{F}=\big{\{}g_{x}(y)=\frac{1}{\sqrt{h^{d}\|K\|_{2}^{2}c}}K\left(\frac{x-y}{h}\right):\;x\in\mathcal{M}\big{\}}, and $\mathbb{G}_{n}(g)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(g(X_{i})-\mathbb{E}g(X_{1})\right),$ $\forall g\in\mathcal{F}.$ we can write $\frac{\sqrt{nh^{d}}(\widehat{f}(x)-\mathbb{E}\widehat{f}(x))}{\sqrt{\|K\|_{2}^{2}c}}=\mathbb{G}_{n}(g_{x}),x\in\mathcal{M}$ and thus

[TABLE]

Let $\mathbb{B}$ be a centered Gaussian process on $\mathcal{F}$ such that for all $g_{x},g_{y}\in\mathcal{F}$, $\mathbb{E}(\mathbb{B}(g_{x})\mathbb{B}(g_{y}))=$ ${\rm Cov}(g_{x}(X_{1}),g_{y}(X_{1})).$ Applying Corollary 2.2 in Chernozhukov et al. (2014) we have that for all $\gamma\in(0,1)$ and $n$ sufficiently large

[TABLE]

where $A_{1},A_{2},A_{3}$ and $A_{4}$ are some constants. See Proposition 3.1 in Chernozhukov et al. (2014) for a similar derivation. Since $\mathbb{E}\left[\sup_{g\in\mathcal{F}}|\mathbb{B}(g)|\right]=O(\sqrt{\log{n}})$ (by Dudley’s inequality for Gaussian processes, c.f. van der Vaart and Wellner, 1996, Corollary 2.2.8), with the choice of $\gamma=1/\log{n}$, we apply Lemma 2.4 in Chernozhukov et al. (2014) and have

[TABLE]

It remains to show that $\lim_{h\rightarrow 0}\mathbb{P}\left(\sup_{g\in\mathcal{F}}|\mathbb{B}(g)|<b(z)\right)=\exp\{-2\exp\{-z\}\}.$ Let $W$ and $B$ be $d$-dimensional Wiener process and Brownian bridge, respectively. Put

[TABLE]

where $M$ is the Rosenblatt transformation (cf. Rosenblatt, 1976). Then

[TABLE]

Let further $\widetilde{U}(x)=\frac{1}{h^{d/2}\|K\|}\int K\left(\frac{x-s}{h}\right)dW(s).$ Following the arguments on page 1013 of Rosenblatt (1976) (also see Proposition 2.2 in Bickel and Rosenblatt, 1973), we have $\sup_{x\in\mathcal{M}}|U(x)-\widetilde{U}(x)|=O_{p}(h^{1/2}).$ We then only need to show

[TABLE]

Next we are going to apply the probability results in Theorem 7.1. Under our assumptions the isosurface $\mathcal{M}$ is a $d-1$ dimensional $C^{1}$ submanifold in $\mathbb{R}^{d}$ (see Theorem 2 in Walther, 1997). As discussed in Remark 2.1b), the reach of the isosurface is positive under our assumptions. It is easy to verify that the conditions for $Q(\delta)$ in (7.3) since $K$ is assumed to have bounded support. We will verify (7.1) and (7.2) in what follows.

First observe that, as $\Delta x/h\rightarrow 0$,

[TABLE]

where $\Sigma$ is a $d\times d$ symmetric matrix with the $(i,j)$-th element

[TABLE]

Note that the little o term in (7.11) is uniform in $x\in f^{-1}[c-\delta_{0},c+\delta_{0}]$ and $h\in(0,1]$, where $\delta_{0}$ appears in assumption (F2). Let $\|\partial K\|_{2}^{2}=\int\left(\frac{\partial K(u)}{\partial u_{1}}\right)^{2}du$. Then due to the symmetry of $K$, $\Sigma_{i,j}=s_{K}^{2}\delta_{i,j}$ where $s_{K}^{2}=\frac{\|\partial K\|_{2}^{2}}{2\|K\|_{2}^{2}}$ and $\delta_{i,j}$ is the Kronecker delta. That is $\Sigma=s_{K}^{2}\bf{I}$. Then (7.1) and (7.2) are satisfied and we can apply Theorem 7.1 to get (7.10) by showing that $b(z)$ is equivalent to

[TABLE]

where $M_{s}$ is a $d\times(d-1)$ matrix with orthonormal columns spanning ${\cal T}_{s}\mathcal{M}$, and $\|\cdot\|_{d-1}$ is the sum of squares of all minors of order $d-1$ for a $d\times(d-1)$ matrix. Note that since $\Sigma=s_{K}^{2}\bf{I}$ we have $\|\Sigma^{1/2}M_{s}\|_{d-1}=s_{K}^{(d-1)}\|M_{s}\|_{d-1}$. By the Cauchy-Binet formula (cf. Broida and Williamson 1989, pp 208-214), $\|M_{s}\|_{d-1}=\sqrt{det(M_{s}^{T}M_{s})}=1$. Recall that $\mathscr{V}_{d-1}$ is the $(d-1)$-dimensional Hausdorff measure. The integral in (7.13) can be simplified as

[TABLE]

Now we have verified the expression in (7.13) is equivalent to $b(z)$ and therefore (2.4) follows.

Next we prove the case $d=1$. Following the discussion after the assumptions, denote $\mathcal{M}=\{x_{i},i=1,\cdots,N\}$. Following similar argument at the beginning of the proof for $d\geq 2$ and using results from Bickel and Rosenblatt (1973), we have that the asymptotic distribution of $\textstyle{\sup\limits_{x\in\mathcal{M}}\frac{\sqrt{nh}|\widehat{f}^{bc}(x)-c|}{\sqrt{\|K\|_{2}^{2}c}}}$ is the same as that of $\sup_{x\in\mathcal{M}}|\tilde{U}(x)|,$ as $n\rightarrow\infty$. When $h$ is small enough, this supremum becomes $\max_{i=1,\cdots,N}Z_{i}$, where $Z_{i}$’s are i.i.d. standard normal random variables. Therefore, as $n\rightarrow\infty$,

[TABLE]

Recall that $\Phi$ is the standard normal c.d.f. Then the c.d.f of $\max_{i=1,\cdots,N}Z_{i}$ is $\Phi^{N}$. By Theorem 3.1 in Biau et al. (2007), $\widehat{N}=N$ for $n$ large enough with probability 1. Following the same argument as for (7.1), we obtain (2.4) for $d=1$. $\square$

7.2 Proof of Corollary 2.1

The proof for $d=1$ is trivial. Now we show the proof for $d\geq 2$. Since

[TABLE]

using Theorem 2.1, we only need to show that

[TABLE]

In fact, we can obtain the following result: for all $L>0$ as $n\rightarrow\infty$,

[TABLE]

Note that the event $\{\widehat{C}_{n,1}^{+}(1-\alpha)\subset\mathcal{L}\subset\widehat{C}_{n,1}^{-}(1-\alpha)\}$ is equivalent to $E_{1}\cap E_{2}$, where $E_{1}=\{f(y)\geq c$, $\forall y$ s.t. $\widehat{f}^{\,bc}(y)\geq c+\widehat{a}_{1-\alpha}^{(d)}$ $\}$, and $E_{2}=\{$ $\widehat{f}^{\,bc}(y)\geq c-\widehat{a}_{1-\alpha}^{(d)}$ $\forall y$ s.t. $f(y)\geq c$ $\}$. We also denote $E_{0}=\{\mathcal{M}\subset\widehat{C}_{n,1}(1-\alpha)\}=\{$ $c-\widehat{a}_{1-\alpha}^{(d)}\leq\widehat{f}^{\,bc}(x)\leq c+\widehat{a}_{1-\alpha}^{(d)}$, $\forall x$ s.t. $f(x)=c$ $\}$. Then (7.16) can be written in the form of $\mathbb{P}(E_{1}\cap E_{2}|E_{0})=1-O(n^{-L})$ or equivalently,

[TABLE]

where $E_{1}^{\complement}=\{$ $\exists y$ s.t. $f(y)<c\;\&$ $\widehat{f}^{\,bc}(y)\geq c+\widehat{a}_{1-\alpha}^{(d)}$ $\}$, and $E_{2}^{\complement}=\{$ $\exists y$ s.t. $f(y)\geq c\;\&$ $\widehat{f}^{\,bc}(y)<c-\widehat{a}_{1-\alpha}^{(d)}$ $\}$. Property (7.17) obviously follows from

[TABLE]

and

[TABLE]

We only show (7.18). To this end, we first introduce two more events. For some $C>0$ large enough, define

[TABLE]

where we use the notation introduced in (2.1). It follows from the proof on page 207 of Mammen and Polonik (2013) that $\mathbb{P}(B_{0})=1-O(n^{-L})$. Also note that

[TABLE]

It is known (e.g. see Theorem 1 in Einmahl and Mason 2005, and Lemma 3 in Arias-Castro et al. 2016 ) that

[TABLE]

Following an argument similar to that in Mammen and Polonik (2013), we thus obtain $\mathbb{P}(B_{1})=1-O(n^{-L})$. These rates of convergences of $P(B_{0})$ and $P(B_{1})$ will be used in the following.

With $\delta_{0}$ given in Assumption (F2), let $\mathcal{M}^{\delta_{0}}=\{x:c-\delta_{0}\leq f(x)\leq c+\delta_{0}\}$. We now split up $E_{1}^{\complement}$ into

[TABLE]

Note that these two sets are disjoint and $E_{1}^{\complement}=E^{\complement}_{11}\cup E^{\complement}_{12}.$ To show (7.18), we now show that both $P(E_{11}^{\complement}|E_{0})=O(n^{-L})$ and $P(E_{12}^{\complement}|E_{0})=O(n^{-L}).$

First we show that $E_{12}^{\complement}\cap E_{0}\cap B_{0}=\emptyset$ for large enough $n$. To this end, let $x\in\mathcal{M}$, i.e. $f(x)=c$. Then, on $E_{0}$, we have $c-\widehat{a}_{1-\alpha}^{(d)}\leq\widehat{f}^{\,bc}(x)\leq c+\widehat{a}_{1-\alpha}^{(d)}$. Let $y$ be a point that makes $E_{12}^{\complement}$ occur, i.e., $y\notin\mathcal{M}^{\delta_{0}}$ and $f(y)<c$ and $\widehat{f}^{\,bc}(y)\geq c+\widehat{a}_{1-\alpha}^{(d)}$. Since $y\notin\mathcal{M}^{\delta_{0}},$ we in fact have $f(y)<c-\delta_{0}$. Note that

[TABLE]

which implies that, on $B_{0}$,

[TABLE]

Also notice that, on $B_{0},$ we have $\widehat{f}^{\,bc}(x)\leq c+C\beta_{n,h}^{(0)}$, and therefore $\widehat{f}^{\,bc}(y)<c-\delta_{0}+3C\beta_{n,h}^{(0)}$, which, for large enough $n$ cannot occur when $\widehat{f}^{\,bc}(y)\geq c+\widehat{a}_{1-\alpha}^{(d)}$. Therefore we get $E_{12}^{\complement}\cap E_{0}\cap B_{0}=\emptyset$ for $n$ large enough.

Next we show that $E_{11}^{\complement}\cap E_{0}\cap B_{1}=\emptyset$ for large enough $n$. So we now assume that $y\in\mathcal{M}^{\delta_{0}}$. Consider the integral curve $\mathfrak{X}_{y}(t)$, which is driven by $\nabla f/\|\nabla f\|^{2}$, starting from $y$. Recall that $\|\nabla f(z)\|>\epsilon_{0}>0$ for $z\in\mathcal{M}^{\delta_{0}}$ by assumption (F2). Let $\theta=c-f(y)>0$ (on $E_{12}^{\complement}$). Using the property of the integral curve described in (3.4), we have $\mathfrak{X}_{y}(\theta)\in\mathcal{M}$. On $B_{1},$ we have $\langle\nabla\widehat{f}^{\,bc}(z),\nabla f(z)/\|\nabla f(z)\|^{2}\rangle>0,\;\forall z\in\{\mathfrak{X}_{y}(t):\;t\in[0,\theta]\}$ for $n$ large enough. This means $\widehat{f}^{\,bc}$ keeps increasing on the trajectory of $\{\mathfrak{X}_{y}(t):\;t\in[0,\theta]\}$. Therefore, $\widehat{f}^{\,bc}(\mathfrak{X}_{y}(\theta))>\widehat{f}^{\,bc}(y)\geq c+\widehat{a}_{1-\alpha}^{(d)}$, which contradicts $E_{0}$. Therefore, $E_{11}^{\complement}\cap E_{0}\cap B_{1}=\emptyset$ for $n$ large enough.

This now results in the following. For $n$ large enough

[TABLE]

Since both $\mathbb{P}(B_{1}^{\complement})=O(n^{-L})$ and $\mathbb{P}(B_{0}^{\complement})=O(n^{-L}),$ and $P(E_{0})\to 1-\alpha$, the assertion follows. The fact that $\mathbb{P}(E_{2}^{\complement}|E_{0})=O(n^{-L})$ can be shown in a similar way. This completes the proof. $\square$

7.3 Proof of Theorem 2.2

Assuming (2.9) (or (2.11)) is true, then (2.10) (or (2.12)) can be proved in a very similar way as for Corollary 2.1. In particular, notice that the key result (7.16) in the proof of Corollary 2.1 can be replaced by (say, for the proof of (2.12))

[TABLE]

Next we prove (2.9) and (2.11). We will use $C$ to denote a generic constant that may be different at different occurrences. For any $\widehat{c}>0$, the following three events are equivalent:

[TABLE]

Similarly

[TABLE]

Therefore it suffices to show that

[TABLE]

and

[TABLE]

The proofs for these two results are similar. We first show (7.21) and then briefly sketch the proof for (7.22). It is known from page 209 in Mammen and Polonik (2013) (also see Theorem 3.1 in Neumann (1998)) that for some $C<\infty$

[TABLE]

for an arbitrarily large $L>0$, which implies

[TABLE]

When $n$ is large enough there exists $C_{1}>0$ such that $\widehat{\mathcal{M}}\subset\bigcup_{x\in\mathcal{M}^{E}}\mathscr{B}(x,C_{1}\beta_{n}^{E})$, where $\mathscr{B}(x,r)$ is the ball in $\mathbb{R}^{d}$ with center $x$ and radius $r$. Therefore

[TABLE]

It follows from an argument similar to the one given in Mammen and Polonik (2013), page 209, that, for some $C>0$,

[TABLE]

This then leads to

[TABLE]

which combining with (7.3) further implies

[TABLE]

where $\gamma_{n}^{E}$ is given in (2.7). As a result of Proposition 3.1 in Neumann (1998) we obtain with $\Psi_{n}$ as in (2.6),

[TABLE]

where this rate holds uniformly in $0\leq c<d<\infty$ Thus, for some $C>0$, we have

[TABLE]

With (7.26) and (7.28), then (7.21) follows from Lemma 2.4 in Mammen and Polonik (2013).

Next, we briefly outline the proof of (7.22). Following the proof on page 207 and Mammen and Polonik (2013) (also see Lemma 3 in Arias-Castro et al., 2016), we have that under our assumption, for some $C>0$,

[TABLE]

Mammen and Polonik (2013, page 210) show that, for all $\delta>0$, there exists $C>0$ such that

[TABLE]

Since $\widehat{\mathcal{M}}\subset\mathscr{B}(0,\delta)$ for large $\delta$ when $n$ is large enough with probability one, by (7.3), (7.29) and (7.3) we have

[TABLE]

Since there exists $C_{2}>0$ such that $\widehat{\mathcal{M}}\subset\bigcup_{x\in\mathcal{M}^{E}}\mathscr{B}(x,C_{2}\beta^{(0)}_{n,h})$ when $n$ is large enough, the following result similar to (7.25) can be derived

[TABLE]

Then combining (7.3) and (7.32) we get

[TABLE]

where $\gamma_{n}$ is given in (2.8). Mammen and Polonik (2013) modify (7.27) to

[TABLE]

which immediately gives, for some $C>0$,

[TABLE]

where (7.34) and (7.35) hold uniformly over $0\leq c<d<\infty$ and $t\in{\mathbb{R}}$, respectively. Applying Lemma 8.1 with $Z_{n}=\sup_{x\in\mathcal{M}}|\widehat{f}(x)-f(x)|$ and $Z_{n}^{*}=\sup_{x\in\widehat{\mathcal{M}}}|\widehat{f}^{*}(x)-\widehat{f}_{g}(x)|$, by using (7.33) and (7.35), we obtain (7.22). $\square$

7.4 Proof of Theorem 2.3

Following the same argument as in the proof of Theorem 2.2, (2.15) follows easily (using the proof of Corollary 2.1) once (2.14) is proved. We will only show the latter. Since

[TABLE]

we have

[TABLE]

Similar to (7.29), there exists $C>0$ such that, for an arbitrary $L>0$,

[TABLE]

and, similar to (7.25), we obtain

[TABLE]

Since \Big{|}\sup_{x\in\widehat{\mathcal{M}}}|\widehat{f}(x)-\mathbb{E}\widehat{f}(x)|-\sup_{x\in\mathcal{M}}|\widehat{f}(x)-\mathbb{E}\widehat{f}(x)|\Big{|}+\sup_{x\in\mathcal{M}}|\widehat{\beta}(x)-\beta(x)|\\ \geq\Big{|}\sup_{x\in\widehat{\mathcal{M}}}|\widehat{f}(x)-\mathbb{E}\widehat{f}(x)|-\sup_{x\in\mathcal{M}}|\widehat{f}^{\,bc}(x)-f(x)|\Big{|}, combining (7.37) and (7.38), we have

[TABLE]

By (7.3), we get

[TABLE]

Similar to (7.34), we have

[TABLE]

where this convergence is uniform in $0\leq c<d<\infty$. Combining (7.37) and (7.41), we have with $\kappa_{n}=Ch^{2}(\beta_{n,l}^{(2)}+h^{2})$,

[TABLE]

This then implies that, for some $C>0$,

[TABLE]

Therefore, using Lemma 8.1 with (7.40) and (7.43), we have

[TABLE]

and the conclusion of the theorem follows. $\square$

7.5 Proofs for Section 3

We begin with two lemmas that will be needed in the proofs. Recall that $\widehat{f}^{\,bc}(x)=\widehat{f}(x)-\widehat{\beta}(x)$ be the de-biased density estimator, and $\widehat{\mathfrak{X}}^{bc}_{x}$ be the integral curve driven by $\frac{\nabla\widehat{f}^{\,bc}}{\|\nabla\widehat{f}^{\,bc}\|^{2}}$ and $\widehat{\theta}^{\,bc}_{x}$ be the corresponding time when $\widehat{\mathfrak{X}}^{bc}_{x}$ hits $\widehat{\mathcal{M}}^{bc}$.

Lemma 7.1

Under assumptions (F1), (F2), (K), $\beta_{n,h}^{(1)}=o(1)$ and $\beta_{n,l}^{(3)}=o(1)$ as $n\rightarrow\infty$, with probability one we have that for $x\in\mathcal{M}$ the solution $\widehat{\theta}_{x}^{bc}$ in (3.6) exists and is unique for $n$ large enough. In such a case the mapping $x\mapsto\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}^{\,bc}_{x})$ is bijective between $\mathcal{M}$ and $\widehat{\mathcal{M}}^{bc}$.

Proof of Lemma 7.1 Throughout this proof we assume that the sample size is large enough. We first show that the solution $\widehat{\theta}_{x}^{bc}$ in (3.6) exists and is unique. Due to the strong consistency of $\widehat{f}^{bc}$, we have that $\widehat{\mathcal{M}}^{bc}\subset\mathcal{M}^{\delta_{0}}:=f^{-1}[c-\delta_{0},c+\delta_{0}]$. Also due to the strong consistency of gradient estimator $\nabla\widehat{f}^{bc}$ implied by (7.19) and (7.20), we have $\|\nabla\widehat{f}^{bc}\|>\epsilon_{0}/2$ for all $x\in\mathcal{M}^{\delta_{0}}$ by assumption (F2). Since the trajectory of $\widehat{\mathfrak{X}}_{x}^{bc}$ is driven by $\nabla\widehat{f}^{bc}$, for $x\in\mathcal{M}$ we have the existence and uniqueness of $\widehat{\theta}_{x}^{bc}$ and $-\infty<\widehat{\theta}_{x}^{bc}<\infty$. Also using $\|\nabla\widehat{f}^{bc}\|>0$ again, as a property of integral curves we have $\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}_{x_{1}}^{bc})\neq\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}_{x_{2}}^{bc})$ if $x_{1}\neq x_{2}$ for $x_{1},x_{2}\in\mathcal{M}$.

The above argument also implies that the mapping $x\mapsto\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}^{\,bc}_{x})$ is bijective between $\mathcal{M}$ and $\widehat{\mathcal{M}}^{bc}$. The injectivity is an immediate consequence. Next we show the surjectivity. For any $y\in\widehat{\mathcal{M}}^{bc}$, without loss of generality we assume $f(y)<c$. Since $\langle\nabla\widehat{f}^{bc},\nabla f/\|\nabla f\|^{2}\rangle>0$ on $\mathcal{M}^{\delta_{0}}$, the value of $f$ keeps increasing on the trajectory of $\widehat{\mathfrak{X}}^{\,bc}_{y}$. Therefore there exists a finite time $\tilde{\theta}_{y}$ when $\widehat{\mathfrak{X}}^{\,bc}_{y}$ hits $\mathcal{M}$. Let $x=\widehat{\mathfrak{X}}^{\,bc}_{y}(\tilde{\theta}_{y})$ and $\widehat{\theta}^{\,bc}_{x}=-\tilde{\theta}_{y}$. Then $y=\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}^{\,bc}_{x})$. $\square$

Lemma 7.2

*Under assumptions (F1), (F2), (K), (H1)2 and (H2)2 we have with

[TABLE]

*that

[TABLE]

and

[TABLE]

Proof of Lemma 7.2 By Lemma 7.1, for $x\in\mathcal{M}$, $\widehat{\theta}_{x}^{bc}$ exists and is unique for large sample. For $x\in\mathcal{M}$, $f(x)=c=\widehat{f}^{\,bc}(\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}^{\,bc}_{x}))$ and therefore

[TABLE]

Consequently,

[TABLE]

where we have used (7.36). This is (7.45). Next we prove (7.46). Without loss of generality we assume $\widehat{\theta}^{\,bc}_{x}>0$ in what follows. We can write

[TABLE]

where we denote $\widehat{\eta}^{\,bc}(t)=\frac{\nabla\widehat{f}^{\,bc}(\widehat{\mathfrak{X}}^{\,bc}_{x}(t))}{\|\nabla\widehat{f}^{\,bc}(\widehat{\mathfrak{X}}^{\,bc}_{x}(t))\|^{2}}-\frac{\nabla\widehat{f}^{\,bc}(\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}^{\,bc}_{x}))}{\|\nabla\widehat{f}^{\,bc}(\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}^{\,bc}_{x}))\|^{2}}$. From (7.5) we obtain

[TABLE]

By (7.47) we have

[TABLE]

Using assumptions (H1)2 and (H2)${}^{2},$ we obtain for some small $\epsilon>0$

[TABLE]

where $\|\cdot\|_{F}$ is the Frobenius norm and $\mathcal{M}\oplus\epsilon=\{x\in\mathbb{R}^{d}:\;d(x,\mathcal{M})\leq\epsilon\}$. By using Taylor expansion, we then have for some small $\epsilon>0$,

[TABLE]

Therefore using (7.5) for some $C>0$, with probability one for $n$ large enough,

[TABLE]

We conclude the proof of (7.46) by using (7.45). $\square$

Proof of Theorem 3.1. First notice that (3.8) (or (3.10)) follows from (3.7 (or (3.9)) using the same approach as in the proof of Corollary 2.1. We point out that the events $E_{1}$ and $E_{2}$ defined in the proof of Corollary 2.1 need to be replaced by (say, for the proof of (3.8))

[TABLE]

Then we can show that

[TABLE]

Details are omitted. We focus on the proofs of (3.10) and (3.9) in what follows.

Part 1. The result immediately follows from Lemma 7.2 and Theorem 2.1.

Part 2. By (7.4) and Lemma 7.2 we obtain that for some $C>0$

[TABLE]

Noting that \tau_{n}^{2}=o\left(\beta_{n,h}^{(1),E}\,\beta_{n,h}^{(0)}+h^{2}\big{(}\beta_{n,l}^{(2)}+h^{2}\,\big{)}\right), the assertion follows from an argument similar to that in the proof of Theorem 2.3. $\square$

The proof of Theorem 3.2. Following the same argument as in the proof of Theorem 3.1, we will only show (3.11) and (3.12) can be derived consequently. Using Lemma 7.2 we have

[TABLE]

Similarly, if we consider the trajectories $\widehat{\mathfrak{X}}_{x}$ traveling between $\mathcal{M}^{E}$ and $\widehat{\mathcal{M}}$, then we have

[TABLE]

Let $\mathbb{P}^{*}$ be the conditional probability measure given $X_{1},\cdots,X_{n}$. Then the bootstrap version of the above result is as follows.

[TABLE]

Notice that $\mathbb{E}\mathbb{P}^{*}=\mathbb{P}.$ Let $\alpha_{n}=C_{1}\left(\beta_{n,g}^{(0)}+h^{2}\right)$ for some $C_{1}$ fixed and large enough. Let $\mathcal{M}\oplus\alpha_{n}=\{x\in\mathbb{R}^{d}:d(x,\mathcal{M})\leq\alpha_{n}\}$. With $\nu_{n,i}$, $i=1,\cdots,5$ to be given in (7.59), (7.60), (7.61), (7.64) and (7.65), we will show the following inequalities:

[TABLE]

Verification of (7.53). It follows from Lemma 7.2 that

[TABLE]

In other words, there exists $C_{0}>0$ such that $\sup_{x\in\mathcal{M}}\|\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}^{\,bc}_{x})-x\|\leq C_{0}\tau_{n}$ for $n$ large enough with probability one. Then we have

[TABLE]

where we have used (7.19) and (7.20). Thus,

[TABLE]

Verification of (7.54). We have that

[TABLE]

Note that

[TABLE]

and

[TABLE]

where the last equality is obtained following a similar argument on pages 208-209 of Mammen and Polonik (2013). Then we have

[TABLE]

Verification of (7.55). Note that

[TABLE]

With (7.3) and (7.37), we have

[TABLE]

Verification of (7.56). First observe that

[TABLE]

We want to show

[TABLE]

On the one hand, $\sqrt{n}[\widehat{f}^{\,*,E}(x)-\mathbb{E}\widehat{f}^{\,*,E}(x)]$ can be viewed as an empirical process indexed by

[TABLE]

Under our assumptions on the bandwidth and kernel function, following a similar argument on page 207 of Mammen and Polonik (2013), we have

[TABLE]

On the other hand, notice that

[TABLE]

A standard argument leads to

[TABLE]

Then (7.62) is a consequence of the above results. Therefore following the same argument as in Cuevas et al. (2006), we have $d_{H}(\widehat{\mathcal{M}}^{\,*,E},\mathcal{M})=O_{a.s.}\left(\alpha_{n}\right),$ which implies that for $C_{1}$ large enough and when the sample size is large enough, $\widehat{\mathcal{M}}^{\,*,E}\subset\mathcal{M}\oplus\alpha_{n}$ with probability one. Therefore

[TABLE]

Following the same argument as for (7.54), we have

[TABLE]

with

[TABLE]

Verification of (7.57). Following an argument similar to (7.58), there exists $C_{0}>0$ such that $\sup_{x\in\widehat{\mathcal{M}}^{*,E}}\|\widehat{\mathfrak{X}}_{x}^{*}(\widehat{\theta}_{x}^{*})-x\|\leq C_{0}\beta_{n,h}^{(0),E}$ when the sample size is large enough with probability one. Then we have

[TABLE]

Here

[TABLE]

Similar to (7.62), we have

[TABLE]

Also, following standard arguments, we obtain

[TABLE]

and

[TABLE]

Therefore

[TABLE]

After collecting the results (7.51), (7.52), (7.53), (7.54), (7.55), (7.56) and (7.57), we have

[TABLE]

because $\gamma_{n}^{bc}+\zeta_{n}/\sqrt{nh^{d}\log{n}}$ is the leading term of \tau_{n}^{2}+\big{(}\beta_{n,h}^{(0),E}\big{)}^{2}+\nu_{n,1}+\nu_{n,2}+\nu_{n,3}+\nu_{n,4}+\nu_{n,5}.

Similar to (7.4), we obtain

[TABLE]

Note that the proof of the above result need to adapt the proof of Proposition 3.1 in Neumann (1998). Specifically, using Neumann’s notation, corresponding to page 2045 in Neumann (1998), we have that

[TABLE]

implies

[TABLE]

where $x_{k}$ is a fixed point on $I_{k}$, say $(2(k_{1}-1)h,\cdots,2(k_{d}-1)h)$, and $\kappa_{n}=Ch\beta_{n,h}^{(0),E}$. This is because $I_{k}$ is a cube with size length of $h$, $\|\nabla f\|$ is differentiable and bounded away from zero in a neighborhood of the level set. Then the rest of the proof follows the proof of Proposition 3.1 in Neumann (1998).

With (7.51), (7.53) and (7.67), we have with $\lambda_{n}=C(\tau_{n}^{2}+\nu_{n,1})$,

[TABLE]

An application of Lemma 8.1 with $Z_{n}=\sup_{x\in\mathcal{M}}\|\widehat{\mathfrak{X}}^{\,bc}_{x}(\widehat{\theta}^{\,bc}_{x})-x\|$ and $Z_{n}^{*}=\sup_{x\in\widehat{\mathcal{M}}^{*,E}}\|\widehat{\mathfrak{X}}_{x}^{*}(\theta_{x}^{*})-x\|$, by using (7.5) and (7.5), concludes the proof. $\square$

8 Appendix

The following result is essentially taken from Mammen and Polonik (2013). We state it here for easy reference.

Using the notation introduced above, let $Z_{n}$ be a statistic, and let $Z_{n}^{*}$ be a bootstrap version of this statistic. For $0<\alpha<1$, define

[TABLE]

and let $c_{n}(1-\alpha)$ be defined similarly with $Z_{n}^{*}$ replaced by $Z_{n}$ (and $P^{*}$ replaced by $P$).

Lemma 8.1

Suppose that there exist sequences $\{\gamma_{n}\},\{\delta_{n}\}$ and $\{\tau_{n}\}$ such that

[TABLE]

Then we have, for $0<\alpha<1,$

[TABLE]

Proof. Lemma 2.4 of Mammen and Polonik (2013) says that under the stated conditions, we have

[TABLE]

It remains to observe that, by using assumption (8.2),

[TABLE]

Note that in Mammen and Polonik (2013), the quantities $Z_{n}$ and $Z_{n}^{*}$ denote particular statistics, their Lemmas 2.2 and 2.4 in fact hold for any statistics satisfying (8.1) and (8.2). Indeed, an inspection of their proofs shows that they do not use any other properties of the statistics. $\square$

Acknowledgement

The authors would like to thank the Associate Editor and the referee for careful reading of the manuscript and for insightful comments that lead to significant improvements. The research of Wolfgang Polonik was partially supported by NSF grant DMS 1107206. The research of Wanli Qiao was partially supported by NSF grant DMS 1821154. The simulations in this work were run on ARGO, a research computing cluster provided by the Office of Research Computing at George Mason University, VA.