Multiscale inference for multivariate deconvolution

In this paper we provide new methodology for inference of the geometric features of a multivariate density in deconvolution. Our approach is based on multiscale tests to detect significant directional derivatives of the unknown density at arbitrary points in arbitrary directions. The multiscale method is used to identify regions of monotonicity and to construct a general procedure for the detection of modes of the multivariate density. Moreover, as an important application a significance test for the presence of a local maximum at a pre-specified point is proposed. The performance of the new methods is investigated from a theoretical point of view and the finite sample properties are illustrated by means of a small simulation study.


Introduction
In many applications such as in biological, medical imaging or signal detection only indirect observations are available for statistical inference, and these problems are called inverse problems in the (statistical) literature. In the case of medical imaging, a well-known example is Positron Emission Tomography. Here, the connection between the 'true' image and the observations involves the Radon transform [see, for example, [10]]. Other typical examples are the reconstruction of biological or astronomical images, where the connection between the true image and the observable image is -at least in a first approximation -given by convolution-type operators [see, for example, [2] or [5]]. Whereas in these models the data is in general described in a regression framework, similar (de-) convolution problems arise in density estimation from indirect observations [see [13] for an early reference]. The corresponding (multivariate) statistical model for density deconvolution is defined by where (Z 1 , ε 1 ), . . . , (Z n , ε n ) ∈ R d × R d are independent identically distributed random variables and the noise terms ε 1 , . . . , ε n are also independent of the random variables Z 1 , . . . , Z n . We assume that the density f ε of the errors ε i is known and are interested in properties of the density f of the random variables Z i based on the sample {Y 1 , . . . , Y n }. In terms of densities, model (1.1) can be rewritten as where g denotes the density of Y 1 . Density estimators can be constructed and investigated similarly to the regression case (see the references in the next paragraph), and in this paper we are interested in describing qualitative features of the density f using the sample {Y 1 , . . . , Y n }. In particular we will develop a method for simultaneous detection of regions of monotonicity of the density f at a controlled level and construct a procedure for the detection of the modes of f . To our best knowledge multivariate problems of this type have not been investigated so far in the literature. On the other hand there exists a wide range of literature concerning statistical inference in the univariate deconvolution model. A Fourier-based estimate of the density f using a damping factor for large frequencies was introduced in [13], whereas [26] estimate f with a wavelet-based deconvolution density estimator [see also [32] for a nonparametric estimator for the corresponding distribution function or [8] for a plug-in estimator of f based on estimation of a scale parameter for the noise level]. [6] develop confidence bands for deconvolution kernel density estimators, while minimax rates for this estimation problem can be found in [9] and [16]. [28] and [18] point out that the detection of regions of monotonicity and of the modes of a density is a more complex problem and [16] shows that the minimax rate for estimating the derivative over a Hölderβ-class (β ≥ 2) in the univariate setting d = 1 is given by n −(β−1)/(2β+2r+1) , where r > 0 denotes the order of polynomial decay of the Fourier transform of the error density f ε . [3] develop a test for the number of modes of a univariate density and [25] proposes a local test for monotonicity for a fixed interval. More recently [30] discuss multiscale tests for qualitative features of a univariate density which provide uniform confidence statements about shape constraints such as local monotonicity properties. These authors use a Komlós-Major-Tusnády (KMT) estimate for the empirical process (cf. [22]). As the classical KMT construction is not suitable for multivariate multiscale problems because it imposes rather strong conditions, it is not obvious how to analyze multiscale inference in a multivariate context. In the present paper we present a solution of this problem. In particular, we use recent results on Gaussian approximations of multivariate empirical processes [ [11]] to address this problem. Multiscale testing is also widely used in spatial testing, see [24] and [31], among others. Here, one aims at the detection of geometric objects of activation in a grid of sensors with noisy measurements and makes use of limit distributions of suprema of sums of i.i.d. Gaussian random variables [cf. e.g. [20]].
Little research has been done regarding multivariate deconvolution problems. Recent references for density estimation are e.g. [12] using kernel density estimators and [29] for a Bayesian approach in the case of an unknown error distribution with replicated proxies available. Hypothesis testing in deconvolution is investigated in [19] and [7].
In the present paper we will develop a multiscale method for simultaneous identification of regions of monotonicity of the multivariate density f in the deconvolution model (1.1). As we do not impose any conditions or even assume prior knowledge about the shape of the density, our problem and approach differ substantially from the methods used in shape-constrained density estimation [see for example [27] and [4], among others, for some references on shapeconstrained density estimation]. In contrast to shape-constrained inference, our approach is based on simultaneous local tests of the directional derivatives of the density f for a significant deviation from zero for "various" directions and locations.
The remaining part of this paper is organized as follows. In Section 2 we present a Fourier based method for the construction of local tests, which will be used for the inference about the monotonicity properties of the density f . Roughly speaking, we propose a multiscale test investigating the sign of the derivatives of the density f in different locations and directions and on different scales. Section 3 is devoted to asymptotic properties, which can be used to obtain a multiscale test for simultaneous confidence statements about the density. Moreover, we also propose a method for the detection and localization of the modes. The finite sample properties of the method are discussed in Section 4 and all proofs are deferred to Sections 5 and 6, while Section 7 contains two technical results.

Multiscale inference in multivariate deconvolution
From a theoretical point of view, only assumptions on the smoothness of φ have to be imposed and therefore, the theoretical part of this paper investigates arbitrary φ. For practical applications, the function φ can be chosen e.g. as a radially symmetric kernel which does not favor any directions, or as a polynomial kernel such as used in the simulations in Section 4. However, this choice has to be made in advance, and φ must be fixed throughout the data analysis. Define For the description of the local monotonicity properties of the function f we introduce the integral If this expression is, say, negative, we can conclude that the derivative of f in direction s has to be strictly larger than zero on a subset of positive Lebesgue Ideally, one would investigate directly the directional derivatives of f for statistical inference regarding its monotonicity properties. However, the estimation of derivatives is difficult, especially in the deconvolution framework of this paper, such that we consider instead the integral (2.1). Note that for h approaching zero the integral (2.1) approximates the directional derivative −∂ s f (t).
In most applications no prior knowledge about the density f is available and therefore, one would like to test for all triples (s, t, h) consisting of all directions s, locations t and scaling factors h. As this is impossible, we choose a finite set of triples T n := {(s j , t j , h j ) | j = 1, . . . , p} and estimate the integral (2.1) for every (s j , t j , h j ) ∈ T n simultaneously. For statistical inference we then propose a multiscale testing procedure. The practical choice of T n depends on the task considered by the experimenter, but typically the choice of an equidistant grid is reasonable. We present below two examples to choose T n to obtain a graphical representation of the density and to obtain a local mode test, respectively.
Statistical inference regarding the monotonicity properties of f is performed by testing simultaneously several hypotheses of the form hj (x) dx, j = 1, . . . , p, see (2.6) below for the estimators. Testing simultaneously means that the same dataset Y i , i = 1, . . . , n, is used for inference about all 2p hypotheses in (2.2) and (2.3), and that we consider the overall error level for at least one false rejection over all tests. To take the multiple testing problem into account, we propose below an investigation of the joint distribution of the p estimates. This approach allows us to control the family wise error rate of the 2p tests for the hypotheses (2.2) and (2.3). Moreover, we can choose p much larger than n, such that standard correction procedures of the p-value in multiple testing problems such as Holm-Bonferroni or False Discovery Rate do not apply.
The method allows for a global understanding of the shape of the density f . A particular feature of the proposed method consists in the fact that by conducting formal statistical tests the multiple level can be controlled (see Theorem 3.2). To be precise, define by T incr n the set of all triples in T n for which the hypothesis (2.2) is rejected, and by T decr n the set of all triples in T n for which the hypothesis (2.3) is rejected. Then the probability of at least one false rejection within the sets T incr n and T decr n can be bounded by a pre-determined error rate α ∈ (0, 1), that is, the method allows to conclude that with probability For example, simultaneous tests for hypotheses of the form (2.2) and (2.3) can be used to obtain a graphical representation of the local monotonicity behavior of the density as displayed in Figure 1 for a bivariate density. at a location t j represents a rejection of the corresponding hypothesis H s j ,t j ,h0 0,incr and provides therefore an indication of a positive directional derivative of f in direction s j at the location t j . For a detailed description of the settings used to provide Figure 1 and an analysis of the results we refer to Section 4.3.
If one is interested in specific shape constraints of the density, say in a test for a mode (local maximum) at a given point x 0 , inference can be conducted investigating the hypotheses for different pairs (t 1 , s 1 ), . . . , (t p , s p ), where t 1 , . . . , t p are points in a neighborhood of x 0 on the lines {x 0 + λs j |λ > 0} (j = 1, . . . , p), respectively (of course, on could additionally use different scales here).
Throughout this paper we will assume that all partial derivatives ∂ s f of the density f are uniformly bounded, such that the estimated quantity (2.1) is bounded by a constant which does not depend on the triple (s, t, h). Using integration by parts, Plancherel's identity and the convolution theorem, we get Here, denote the Fourier transform and its inverse, respectively, z is the complex conjugate of z ∈ C and x.y stands for the standard inner product of x, y ∈ R d .
For the construction of tests for the hypotheses in (2.2) and (2.3) we define the statistic Because (by (2.5))

Asymptotic properties
In this section we investigate the asymptotic properties of a statistic which can be used to control the multiple level of the tests introduced in Section 2. To be precise, we consider the finite subset of cardinality p ≤ n K for the calculation of the maximum of appropriately standardized statistics T n s,t,h , where K > 1 and for some ε > 0 Throughout this paper we will make frequent use of multi-index notation, where α = (α 1 , . . . , α d ) ∈ N d 0 denotes a multi-index (written in bold), |α| = α 1 + . . . + α d its "length", and for a sufficiently smooth function f : R d → R and a multiindex α we denote by its partial derivative.

K. Eckle et al.
Recall the definition of F s,t,h in (2.7) and to simplify the notation define for a point ( The testing procedure for the hypotheses (2.2) and (2.3) is based on the p estimates T n s j ,t j ,hj = 1 For a rigorous statistical test which controls the multiple level we therefore need to investigate the asymptotic joint distribution of Recall that p is growing with n. Thus, the maximum over all p random variables in (3.3) is in general not bounded. As a consequence, the random variables defined in (3.3) have to be properly standardized. It turns out that the approriate standardization is given bỹ (3.4) whereĝ n is a density estimator of g satisfying g −ĝ n ∞ = o(log(n) −1 ) almost surely (3.5) (for example a kernel density estimator as considered in [17]) and The quantity V j is well-defined under the assumptions presented below (see Lemma 5.2 for details).
Note that the boundary of the hypotheses H s j ,t j ,hj 0,incr and H s j ,t j ,hj 0,decr in (2.2) and hj (x) dx = 0 and in this case we have Consequently, we will investigate the asymptotic properties of max 1≤j≤pX (1) j in the following discussion. For this purpose we make the following assumptions. Assumption 1. Assume that the density g is Lipschitz continuous and locally bounded from below, i.e.

Assumption 2.
We assume a polynomial decay of the Fourier transform of the error density f ε , i.e. that there exist constants r > 0 for d ≥ 2 resp. r > 1 2 for d = 1 and 0 < C u < C o such that Note that as a direct consequence of Assumption 1 g is bounded from above and that there exists a constant δ > 0 such that g(x) ≥ c 2 > 0 for all x ∈ [−δ, 1 + δ] d . Assumption 2 can be seen as a multivariate generalization of the classical assumptions on the decay of the Fourier transform of the error density in the ordinary smooth case (see e.g. [30], Assumption 2). We also note that this assumption defines a mildly ill-posed situation (see [7]). The next assumptions refer to the kernel φ and are required for some technical arguments.

Assumption 3.
Let ∂ s φ L 2 (R d ) = 0 for all s ∈ S d−1 and assume that ∂ β φ exists in [−1, 1] d and is continuous for all |β| ≤ r + 2 , where r is the constant from Assumption 2. We assume further that for some δ > 0 the inequality As for all s ∈ S d−1 and some constant C > 0 that only depends on d, Assumption 3 yields a uniform upper bound for the integral Recall the definition ofX (1) j in (3.4) and define the vector Our first main result provides a uniform approximation of the probabilities P(X (1) ∈ A) by the probabilities P(X ∈ A) for every half-open hyperrectangle A, where the components of the vectorX = (X 1 , . . . ,X p ) are defined by The limit processX is constructed in such a way that it has (asymptotically) the same covariance structure as the vectorX (1) consisting of all test statistics. Moreover, the processX j does not depend on unknown quantities. In order to construct quantiles for the testing procedure for the hypotheses (2.2) and (2.3), we consider the quantity max 1≤j≤pXj and use Theorem (3.1) below to show that the quantiles of max 1≤j≤pX (1) j can be approximated by those of max 1≤j≤pXj (note that this is a simple consequence of Theorem 3.1, using the

Theorem 3.1. Let A denote the set
Then, Furthermore, the random variable max 1≤j≤pXj is almost surely bounded uniformly with respect to n.
Theorem 3.1 will be used to control the multiple level of statistical tests for the hypotheses of the form (2.2) and (2.3). To this end, let α ∈ (0, 1) and denote by κ n (α) the smallest number such that (3.9) By Theorem 3.1, κ n (α) is bounded uniformly with respect to n. The jth hypothesis in (2.2) is rejected, whenever Note that the ill-posedness of the deconvolution problem is reflected in the value of the quantile κ j n (α) through the multiplication with h −r j and by the standardization with the quantity V j .

Theorem 3.2.
Assume that the tests (3.10) and (3.12) for the hypotheses (2.2) and (2.3) are performed simultaneously for j = 1, . . . , p. The probability of at least one false rejection of any of the tests is asymptotically at most α, that is for n → ∞.
It is a well-known fact that statistical inference regarding the qualitative features of a multivariate density is a challenging problem from a computational point of view. In the present context conducting all tests (3.10) and (3.12) for the hypotheses (2.2) and (2.3) is computationally demanding. In general, the support of the deconvolution kernel F j is not compact and therefore, the computation of all p test statistics consists of p · n kernel evaluations.
..,p of the Gaussian limit process depends on p · (p + 1)/2 numerical integrations and for the determination of the quantiles of the limit process p-dimensional normal distributed random vectors have to be simulated.
Next we introduce a method for the detection and localization of the modes of the density. The main idea is to conduct the local tests for modality proposed in (2.4) for a set of candidate modes which does not assume any prior knowledge about the density. To be precise, we assume the following condition on the set T n : for any fixed h and s the set {t : (s, t, h) ∈ T n } is an equidistant grid in [0, 1] d with grid width h. Furthermore, for any fixed t and h the set {s : (s, t, h) ∈ T n } is a grid in S d−1 with grid width converging to zero with increasing sample size.
This grid is now used as follows to check if a point dh for some c > 2 √ d sufficiently large and angle(t − x 0 , s) → 0 for n → ∞. By the condition on T n defined above, the set T x 0 n is nonempty for sufficiently large n. We now use the local tests (3.12) for the hypotheses (2.4) and decide for a mode at the point x 0 if the null hypotheses in (2.4) are rejected for all triples in T x 0 n simultaneously. Note that by choosing the test locations as the vertices of an equidistant grid no prior knowledge about the location of x 0 has to be assumed. Theorem 3.3 below states that the procedure detects all modes of the density with asymptotic probability one as n → ∞. Theorem 3.3. Let x 0 ∈ (0, 1) d denote an arbitrary mode of the density f and assume that there exist functions If the set for some C > 0 sufficiently large is nonempty, then the procedure described in the previous paragraph detects the mode x 0 with asymptotic probability one as n → ∞.
The method to detect the modes of the density proposed in Theorem 3.3 proceeds in two steps: the verification of the presence of a mode with asymptotic probability one in the asymptotic regime presented above and its localization at the rate n −1/(d+2r+4) (up to some logarithmic factor) given by the grid width. [16] showed that in the univariate setting d = 1 the minimax rate for estimating the derivative of a density in a deconvolution problem over a Hölder-β-class is of order n −(β−1)/(2β+2r+1) (β ≥ 2), and it is conjectured that the rate is of order n −(β−1)/(2β+2r+d) in the multivariate case. In the case of mode estimation there are no results available regarding optimal rates of estimates (to the best of our knowledge). However, as the problem of estimating a derivative is closely related to mode estimation, we expect similar optimal rates in the context considered in this paper. In the case β = 2 the optimal rate for estimating the derivative is n −1/(d+2r+4) and Theorem 3.3 shows that the proposed mode estimator attains this rate up to a logarithmic factor. An important and challenging problem for future research is to prove that these rates are in fact minimax optimal.

Finite sample properties
In this section we illustrate the finite sample properties of the proposed multiscale inference. The performance of the test for modality at a given point x 0 (see the hypotheses in (2.4)) and the dependence of its power on the bandwidth and the error variance is investigated. We also illustrate how simultaneous tests for hypotheses of the form (2.2) and (2.3) can be used to obtain a graphical representation of the local monotonicity properties of the density.
We consider two-dimensional densities, i.e. d = 2. The density f ε of the errors in model (1.1) is given by a symmetric bivariate Laplacian with scale parameter σ > 0 which is defined through its characteristic function for (y 1 , y 2 ) ∈ R 2 (cf. [23], Chapter 5). This means that r = 2 and straightforward calculations show that The test function is chosen as where c 2 defines the normalization constant, that is (note that φ is smooth within its support). Moreover, the integration by parts formula gives as φ vanishes on the boundary of its support. Finally, by the representation (4.2) we find that the deconvolution kernel possesses all properties that are used for the proof of Theorem 3.1 and therefore Theorem 3.1 is also satisfied for the function φ. Throughout this section the nominal level is fixed as α = 0.05, and level and power are always stated in percent.

A local test for modality -testing for a single mode
In this section we investigate the performance of a local test for the existence of a mode (more precisely a local maximum) at a given location x 0 which is defined by testing several hypotheses of the form (2.4) simultaneously. Moreover, the influence of the choice of the different parameters on the power of the test is also investigated. To be precise, we conduct four tests for the hypotheses (2.4) with a fixed bandwidth h = h 0 . The postulated mode is given by the point x 0 = (0, 0) and the four directions and locations are chosen as s 1 = t 1 = (1, 0) , s 2 = t 2 = (0, 1) , s 3 = t 3 = (−1, 0) and s 4 = t 4 = (0, −1) . We conclude that f has a local maximum at the point x 0 = (0, 0) , whenever all hypotheses where κ j n (α) is defined by (3.11). An illustration of the considered situation is provided in Figure 2. The quantiles κ n (0.05) defined in (3.9) are derived by 1000 simulation runs based on normal distributed random vectors. In Table  1 we display the normalized quantiles √ nκ 1 n (0.05) for the sample sizes n = 500, 1000, 4000 observations and h 0 = 0.5. Here, the value of the parameter of the Laplacian error density has been chosen as σ = 0.075.
The approximation of the level of the test for a mode at the point x 0 defined by (4.3) is investigated using a uniform distribution on the square [−2.5, 2.5] 2 for the density f . Recall that the quantiles κ j n (α) are constructed in such a way that the probability of at least one false rejection of any of the tests (4.3) is at most α. However, the detection of the presence of a mode is based on simultaneous   The simulated rejection probabilities are presented in Table 3 and show that the mode test performs well, even for small sample sizes. We further note the superiority of the calibrated test. Moreover, we find that the shape of the modal region, which is determined by the absolute values of the eigenvalues of the covariance matrix, has a strong influence on the power of the test (4.3). In the case of N (0, Σ 1 )-distributed random variables Z i (eigenvalues approximately 1.8 and 0.3) the test performs better as for standard normal observations (with both eigenvalues equal to one). In the case of N (0, Σ 2 )-distributed random variables Z i (eigenvalues approximately 3.4 and 0.3) the test performs slightly worse than in the first case but still better as for standard normal observations due to the eigenvalue with absolute value smaller than one. Dependence of the power on a misspecification of the position of the mode: We also investigate the influence of a (slight) misspecification of the position of the candidate mode on the power of the test (4.3) in the situation considered in Table 3 with covariance matrix I for the candidate mode x 0 = (0.2, 0.2) . The results are presented in Table 4 and should be compared with the second and third column in Table 3. We find that the slight misspecification of the position of the candidate mode affects the power of the method only slightly. Dependence of the power on the bandwidth: Next we fix the number of observations, that is n = 1000, the value of the parameter σ = 0.075 and vary the bandwidth h 0 to investigate its influence on the power of the test (4.3). Recall that by the proposed choice of a Laplacian error density, the deconvolution kernel has compact support in [−1, 1] 2 . Hence, by dividing the bandwidth by 2 a fourth of the area is considered and (roughly) a fourth of the number of observations is used for the local test. Thus, we observe a decrease in power of the test for decreasing values of bandwidths which is illustrated in Table 5 Table 6 and we observe that an increase in the value of σ decreases the power of the test. On the other hand the power of the test is very stable for small values of σ.  Table 6 Dependence of the power of the test (4.3) for a mode at the point x 0 = (0, 0) on the scale parameter in the situation considered in Table 3

A local test for modality -testing for two modes simultaneously
We also consider a bimodal density and conduct simultaneously local tests for modality based on the hypotheses (2.4) for the candidate modes x 1 = (0, 0) and x 2 = (3, 0) . We conduct eight tests for the hypotheses   where the quantile κ j n (α) is defined by (3.11). An illustration of the considered scales is provided in Figure 3.  Table 7.  Table 7 Simulated level and power of the tests (4.5) and We observe that in the symmetric case the test detects both modes with (roughly) the same power, whereas in the asymmetric case the mode with smaller variance (even though there is a slight misspecification of its position) is detected more often.
A scatter plot of n = 4000 observations from the convolution of the asymmetric bimodal density and a bivariate Laplace distribution with scale parameter

Inference about local monotonicity of a multivariate density
The multiscale approach introduced in Section 2 can be used to obtain a graphical representation of the monotonicity behavior of a (bivariate) density. We construct a global map indicating monotonicity properties of the density f by conducting the tests (3.10) for the hypotheses (2.2) for a fixed bandwidth of h = 0.5. The set of test locations T t is defined as the set of vertices of an equidistant grid in the square [−1, 2] 2 with width 1 and the set of test directions is given by The tests (3.10) are conducted for every triple The scaling factor for the Laplace density in the convolution model (1.1) is given by σ = 0.075. We consider the tri-modal density with differently shaped modal regions displayed in Figure 5. Figure 1 in Section 2 provides the graphical representation of the monotonicity behavior of the density f . Here, each arrow at a location t in direction s displays a rejection of a hypothesis (2.2). The map indicates the existence of modes close to the points (−0.5, −0.5) , (1.5, −0.5) and (0.5, 1.5) .

Proof of Theorem 3.1
We split the proof of Theorem 3.1 in three parts. The first part is dedicated to several auxiliary results involving the deconvolution kernel F s,t,h . In the second part of the proof we show the approximation (3.8). Finally we conclude by proving the boundedness of the limit distribution in the third part.
Throughout this section the symbols and mean less or equal and greater or equal, respectively, up to a multiplicative constant independent of n and (s, t, h), and the symbol |a s,t,h | |b s,t,h | means that |a s,t,h /b s,t,h | is bounded from above and below by positive constants.

Auxiliary results
We begin with some basic transformations of the deconvolution kernel F s,t,h . Recall that by definition of the kernel φ t,h and the Fourier transform. A substitution in the inner integral shows that By the definition of the inverse Fourier transform and a substitution in the outer integral, we obtain where i denotes the imaginary unit. The following lemma presents some immediate consequences of the Assumptions 2 and 3 made in Section 3.
Proof of Lemma 5.1. 1.: An application of Cauchy-Schwartz's inequality yields for any δ > 0 By Assumption 3, there exists a constant δ > 0 such that the latter integral is bounded uniformly with respect to s. Hence, the assertion follows from the integrability of the function (1 + y 2 ) −(d+δ)/2 .

2.: By Leibniz's rule we have
Moreover, from Lemma 7.2 it follows that where M k is the set of all k-tuples of non-negative integers satisfying k j=1 jm j = k. Assumption 2 in Section 3 yields the estimates Thus, as k j=1 jm j = k for some (m 1 , . . . , m k ) ∈ M k , we find In the case r ≥ k, the claim is now a direct consequence of the estimate similar arguments as given in proof of 1. and Assumption 3.
If r < k we divide the integration area into the ball B 1 (0) and its complement. For the integral h −r . Therefore, we can bound the integral over the complement of the unit ball by the integral over R d and proceed similarly to the first case. It remains to consider the integral over the ball B 1 (0). To this end, notice that Hence, by the boundedness of ∂m −k ∂ym −k l F (∂ s φ) (which follows from the compactness of the support of φ) it remains to show that the integral is bounded, where we used a polar coordinate transform to obtain the inequality. As k ≤ (d + 1)/2 and r > 0, the integral on the right hand side is obviously finite.
Part 1 of the following lemma shows that the constants V 1 , . . . , V p defined in (3.6) are uniformly bounded from above and below.

Lemma 5.2. It holds
3) It now follows from Assumption 2 and a substitution that and the latter integral is bounded by Assumption 3 which concludes the proof of the upper bound.
For the lower bound we find from (5.3) and Assumption 2 that for any constant a > 0. Moreover, for a sufficiently small radius a by the integrability of |F (∂ s φ)| 2 (Assumption 3) and Plancherel's theorem. Furthermore, the mapping s → ∂ s φ L 2 (R d ) is continuous such that by Assumption 3 ∂ s φ L 2 (R d ) ≥ c > 0 for a constant c that does not depend on s.

2.: The representation (5.2) and a substitution in the integral for the variable
As

the differentiation rule for Fourier transforms yields
where the last identity follows from Plancherel's theorem. We now proceed similarly as in the proof of Lemma 5.1 2. and note that An application of the Assumptions 2 and 3 shows Moreover, by Assumption 2, we have This concludes the proof for r ≥ 1. For r < 1 we split up the area of integration into the ball B 1 (0) and its complement and find the required result for the integration over the complement using similar arguments as in the proof of Lemma 5.1 2. For the integral over the unit ball we also follow the line of arguments presented in the proof of Lemma 5.1 2. which yields the required result provided that the integral on the right hand side of the inequality exists. This is the case for all r > 0 if d ≥ 2 and all r > 1 2 in the case d = 1. 3. and 4.: These are direct consequences of Hölder's inequality and 1. resp.

2.
The following Lemma will be used in the second part of the proof of Theorem 3.1.

Proof of Lemma 5.3. 1.: Using the representation (5.2) and Assumption 2 it follows that
The claim follows from the uniform boundedness of S s j shown in Lemma 5.1 1.
2.: Using the representation (5.2), the boundedness of the density g and a substitution we get The proof will be completed showing the estimate For this purpose we decompose the domain of integration for the variable x in two parts: the cube [−δ, δ] d for some δ > 0 and its complement. For the integral with respect to the cube we use the upper bound R d provided in the proof of 1. which yields the required result. For the integral with respect to ([−δ, δ] d ) C note that where the sets A k,l are defined by

Nowm = (d + 1)/m fold integration by parts yields
as |x l | ≤ |x l | for all l = l and |x l | > δ in A k,l .

Proof of the approximation (3.8)
For the consideration of the absolute values we introduce the set and denote by A the set of all hyperrectangles in R 2p of the form We will show below in Section 5.2.1 that the random vectors for any q > 0, where Y 1 , . . . , Y n are independent random vectors, Note that we have , as the random variables X 1 , . . . , X n are i.i.d. and Y 1 , . . . , Y n are independent. Introduce a Gaussian process (B(Φ)) Φ∈L ∞ (R d ) indexed by L ∞ (R d ) as a process whose mean and covariance functions are 0 and respectively. Hence, there exists a version ofB(Φ) such that To derive an alternative representation of the processB recall the definition of the isonormal process (B(Φ)) Φ∈L 2 (R d ) as a Gaussian process whose mean and covariance functions are 0 and R d Φ 1 (x)Φ 2 (x) dx, respectively (see, e.g. [21], Section 5.1). In particular, note that (B(1 A )) A∈B(R d ) defines white noise, where B(R d ) denotes the Borel-σ-field on R d . Throughout this paper, we will use the notation There exists a version of the isonormal process such thatB defines a Gaussian process with the covariance kernel (5.5)). Thus, .
uniformly with respect to s, t, h (by assumption). Furthermore, max .
An application of Markov's inequality finally proves Here, we have investigated convergence in probability w.r.t. the sup-norm. However, standard arguments show that this implies the convergence which is investigated in Theorem 3.1.
In a second step we find that the normalization with c j := ( g(t j )V j ) −1 , j = 1, . . . , 2p, has no influence on the convergence as translation and multiplication preserve the interval structure. More precisely, for any set still defines an element of the set A . A similar result holds for the normalization of the test statistic.
In a third step we show in Section 5.2.2 that the normalization with the density estimator yields to a distribution-free limit process. We firstly assume that the density g is known and prove Hence, by the consideration of the symmetric set T n it follows from (5.4), (5.7) and (5.9) that as for any real valued random variable X and any a ∈ R it holds Next we insert the bandwidth normalization terms. To this end, we introduce the notation and write w j = w(h j ),w j =w(h j ). Similar arguments as in (5.8) show that the insertion of the bandwidth correction terms has no influence on the convergence. Thus recalling the definition ofX j = w j h 11) and it remains to replace the true density by its estimator. For this purpose we show that almost surely by the boundedness from below of g (and therefore ofĝ n almost surely). A null addition of the termw j shows that the latter is equal to The claim follows now from the convergence of w j proven in (5.11) and the a.s. boundedness of the maximum of the limiting process proven in Section 5.3 below. Note that we used the fact that is decreasing in a neighborhood of 0 (cf. [30], Lemma B.11).

Proof of (5.4)
The proof of (5.4) mainly relies on Proposition 2.1 in [11]. The result is stated as follows.
. . , n. Let b, q > 0 be some constants and let B n ≥ 1 be a sequence of constants, possibly growing to infinity as n → ∞. Assume that the following conditions are satisfied: Then, where the sequences D (1) n and D (2) n,q are given by , D (2) n,q = For an application of Theorem 5.4 we have to verify the condition 1. and to find an appropriate sequence B n for conditions 2. and 3. For a proof of condition 1. notice that where we used (5.6) in the inequality. Moreover, as the density of g is bounded from below (Assumption 1) we have In Lemma 5.2 1. we have proven that F j , and using the representation (5.2) we obtain We now follow the line of arguments presented in the proof of Lemma 5.3 2. for m = 2 and note that by conducting integration by parts we get an additional factor h d+1 j . Hence, This concludes the proof of condition 1. as E(X 2 For a proof of condition 2. note that by part 1 of Lemma 5.3 it follows that for any q > 0, which proves (5.4).

5.2.2.
Proof of (5.9) Define then the assertion follows from the statement Here, we used the fact that the constants V 1 , . . . , V 2p are bounded uniformly from below (cf. Lemma 5.2). For this purpose, we will make use of a Slepiantype result. Note that for all δ > 0 (5.14) For the first integral on the right hand side of (5.14) we use the Lipschitz continuity of g (Assumption 1) and find for some ξ satisfying |ξ −g(t j )| ≤ |g(x)−g(t j )|. If δ > 0 is sufficiently small, then g is bounded from below on [−δ, 1 + δ] d (see the remark following Assumption 1), and Lemma 5.2 2. shows that an upper bound of this term (up to some constant) is given by The second integral on the right hand side of (5.14) is bounded by h max which follows from (5.12) and the boundedness of g (Assumption 1). Summarizing, we obtain E(R 2 j ) h max . Moreover, we can show by similar calculations as presented above and an application of Lemma 5.2 4. that Introducing the random variables we obtain from Lemma 5.2 1. and 3.
Hence, max and Theorem 2.2.5 in [1] yields Note that by the symmetry of the set T n with respect to the direction we have and we can consider expectations of positive random variables here.
For an upper bound of E(max 1≤j≤2pRj ) we use the a.s. asymptotic boundedness of shown in Section 5.3 below, which implies and therefore E(max 1≤j≤2p R j ) = O( h max log(n)). This proves (5.9) by an application of Markov's inequality.

Boundedness of the approximating statistic
In order to prove that the approximating statistic max 1≤j≤pXj considered in Theorem 3.1 is almost surely bounded uniformly with respect to n ∈ N we note that for all p ∈ N max where the random variable B is defined by . B does not depend on n and we show below that B is almost surely bounded. We will make use of the following result (Theorem 6.1 and Remark 1, [14]).
Theorem 5.5. Let X be a stochastic process on a pseudometric space (T , ρ) with continuous sample paths. Suppose that the following three conditions are satisfied.
2. For some constants L, M ≥ 1, is finite almost surely.

For some constants
For the application of Theorem 5.5 we introduce the pseudometric space In the following, we prove that the process X fulfills the conditions of Theorem 5.5.

1.:
We have by definition of σ and ρ that as X(s, t, h)/σ(h) corresponds in distribution to a normal distributed random variable with mean zero and variance one by definition of V s,t,h .
2.: By definition, X(s 1 , t 1 , h 1 ) − X(s 2 , t 2 , h 2 ) corresponds in distribution to a normal distributed random variable with mean zero and variance h d+r+1 .
W.l.o.g. we assume in the following h 1 ≤ h 2 and note that condition 2. (with L = 2) follows from the inequality In the first inequality we used the fact that V s 1 ,t 1 ,h1 is uniformly bounded from below and In a proof of the second inequality in (5.15) we note that by application of the triangle inequality Moreover, we find by another application of the inequality F s 2 ,t 2 ,h2 Hence, observing (5.16) and (5.17) the inequality (5.15) follows from ). (5.18) For a proof of this inequality we use Plancherel's theorem which yields The integrand on the right hand side can be estimated as follows and we obtain F s 1 ,t 1 ,h1 − F s 2 ,t 2 ,h2 where e k denotes the kth unit vector of R d (k = 1, . . . , d). By a substitution it follows that Here, we used another substitution and the triangle inequality. For an upper bound for the first term on the right hand side of (5.19), note that by Assumption 3 R d (1+ y 2 ) r |F (∂ e k φ)(y)| 2 dy is finite. Furthermore, a substitution within the Fourier transform shows that the second term of the right hand side of (5.19) is not greater than By an application of Euler's formula, cos(x) ≥ 1 − x for all x ≥ 0 and Cauchy-Schwartz's inequality, we find Therefore, two substitutions and Assumption 3 show that the second term on the right hand side of (5.19) is bounded from above (up to some constant) by It remains to consider the third term on the right hand side of (5.19). Plancherel's theorem, the rule for the Fourier transform of a derivative and a substitution show that the third term on the right hand side of (5.19) can be bounded by |α|≤ r+1 Observe that the distance between adjacent points in the set T 3 := (jε) 1/d , j = 1, . . . , δ ε is equal to ε. As a consequence,Ñ (ε, (0, δ 1/d ]) δ ε . From (5.21) and the results presented above we deduce N (δu) where the latter inequality has been proven in 2. Hence, similar arguments as presented in 3. show thatÑ (ε, T ,d) ε −a for some a > 0, which concludes the proof of the a.s. continuity of the sample paths of X w.r.t.d and implies the a.s. continuity of the sample paths of X w.r.t. ρ.

Proof of Theorem 3.2.
Denote by q the probability of at least one false rejection among all tests (3.10) and (3.12). Using Theorem 3.1, we further deduce from (3.9) for n → ∞.

Proof of Theorem 3.3.
We begin deriving a criterion for the simultaneous rejection of the hypotheses (2.3) on a given set of scales. To this end, let 0 < (α n ) n∈N < 1 be an arbitrary null sequence and J ⊆ {1, . . . , p} be the set of all indices where the inequality hj (x) dx > 2κ j n (α n ) (6.1) is satisfied. An application of Theorem 3.1 shows that the probability of simultaneous rejection of the Null Hypotheses for all tests in (3.12) indexed by J (where α is replaced by α n ) is asymptotically equal to one, i.e.
q := P n −1 by similar arguments as presented in the proof of Theorem 3.2. Now let x 0 ∈ (0, 1) d be a mode of f and (s, t, h) ∈ T x 0 n , i.e. ch ≥ x 0 − t ≥ 2 √ dh for some c > 2 √ d and angle(x 0 − t, s) → 0 for n → ∞. Following the line of arguments presented in the proof of Theorem 3.3 in [15], one can prove that, under the given assumptions, ∂ s f (x) −h for all x ∈ suppφ t,h . Hence, By Theorem 3.1, we find that For a proof of (6.1) it remains to find a condition on h such that which holds for h ≥ C log(n) 1/(d+2r+4) n −1/(d+2r+4) for some C > 0 sufficiently large.