Simultaneous Confidence Bands for Functional Data Using the Gaussian Kinematic Formula

TL;DR

This paper develops a method for constructing simultaneous confidence bands for functional data using the Gaussian Kinematic formula, demonstrating its effectiveness even with small samples and extending it to noisy, discretely sampled data.

Contribution

It introduces a CLT-based approach for applying the Gaussian Kinematic formula to non-Gaussian functional data, with practical extensions and a comprehensive simulation study.

Findings

tGKF performs well with small sample sizes

tGKF outperforms bootstrap methods for larger samples

Method is computationally efficient and applicable to DTI data

Abstract

This article constructs simultaneous confidence bands (SCBs) for functional parameters using the Gaussian Kinematic formula of $t$-processes (tGKF). Although the tGKF relies on Gaussianity, we show that a central limit theorem (CLT) for the parameter of interest is enough to obtain asymptotically precise covering rates even for non-Gaussian processes. As a proof of concept we study the functional signal-plus-noise model and derive a CLT for an estimator of the Lipschitz-Killing curvatures, the only data dependent quantities in the tGKF SCBs. Extensions to discrete sampling with additive observation noise are discussed using scale space ideas from regression analysis. Here we provide sufficient conditions on the processes and kernels to obtain convergence of the functional scale space surface. The theoretical work is accompanied by a simulation study comparing different methods to construct SCBs for the population mean. We show that the tGKF works well even for small sample sizes and only a Rademacher multiplier-$t$ bootstrap performs similarily well. For larger sample sizes the tGKF often outperforms the bootstrap methods and is computational faster. We apply the method to diffusion tensor imaging (DTI) fibers using a scale space approach for the difference of population means. R code is available in our Rpackage SCBfda.

Figures

Tables

Table 1. Table 1 : Average of estimates q ^ α , N subscript ^ 𝑞 𝛼 𝑁 \hat{q}_{\alpha,N} of the width and twice its standard error for different methods of SCBs and the Gaussian Model B 𝐵 B . The simulations are based on 1,000 Monte Carlo simulations.

N	$\mathrm{𝟏𝟎}$	$\mathrm{𝟐𝟎}$	$\mathrm{𝟑𝟎}$	$\mathrm{𝟓𝟎}$	$\mathrm{𝟏𝟎𝟎}$	$\mathrm{𝟏𝟓𝟎}$
true	4.118	3.382	3.211	3.081	2.993	2.96
tGKF	$3.980\pm 0.033$	$3.368\pm 0.011$	$3.204\pm 0.006$	$3.084\pm 0.004$	$3.000\pm 0.002$	$2.973\pm 0.001$
GKF	$2.929\pm 0.015$	$2.925\pm 0.008$	$2.923\pm 0.005$	$2.922\pm 0.003$	$2.921\pm 0.002$	$2.921\pm 0.001$
Degras	$2.866\pm 0.021$	$2.881\pm 0.010$	$2.887\pm 0.007$	$2.890\pm 0.005$	$2.892\pm 0.003$	$2.896\pm 0.003$
Boots	$2.690\pm 0.022$	$2.813\pm 0.013$	$2.851\pm 0.010$	$2.878\pm 0.007$	$2.897\pm 0.005$	$2.903\pm 0.004$
Boots-t	$5.728\pm 0.507$	$3.609\pm 0.050$	$3.292\pm 0.020$	$3.109\pm 0.009$	$2.998\pm 0.005$	$2.968\pm 0.004$
gMult	$3.020\pm 0.161$	$2.921\pm 0.032$	$2.914\pm 0.017$	$2.910\pm 0.008$	$2.913\pm 0.005$	$2.913\pm 0.004$
gMult-t	$3.428\pm 0.046$	$3.091\pm 0.017$	$3.011\pm 0.011$	$2.960\pm 0.007$	$2.931\pm 0.004$	$2.925\pm 0.004$
rMult	$2.690\pm 0.179$	$2.766\pm 0.032$	$2.808\pm 0.015$	$2.849\pm 0.009$	$2.882\pm 0.005$	$2.893\pm 0.004$
rMult-t	$4.100\pm 0.122$	$3.382\pm 0.023$	$3.202\pm 0.013$	$3.078\pm 0.008$	$2.993\pm 0.005$	$2.965\pm 0.004$

Table 2. Table 2 : Average of estimates q ^ α , N subscript ^ 𝑞 𝛼 𝑁 \hat{q}_{\alpha,N} of the width and twice its standard error for different methods of SCBs and the non-Gaussian Model B 𝐵 B with ν = 7 𝜈 7 \nu=7 . The simulations are based on 1,000 Monte Carlo simulations.

N	$\mathrm{𝟏𝟎}$	$\mathrm{𝟐𝟎}$	$\mathrm{𝟑𝟎}$	$\mathrm{𝟓𝟎}$	$\mathrm{𝟏𝟎𝟎}$	$\mathrm{𝟏𝟓𝟎}$
true	4.532	3.628	3.373	3.19	3.048	3.004
tGKF	$3.973\pm 0.040$	$3.365\pm 0.013$	$3.203\pm 0.008$	$3.083\pm 0.004$	$3.000\pm 0.002$	$2.973\pm 0.001$
GKF	$2.926\pm 0.019$	$2.922\pm 0.009$	$2.922\pm 0.006$	$2.921\pm 0.004$	$2.921\pm 0.002$	$2.921\pm 0.001$
Degras	$2.864\pm 0.021$	$2.880\pm 0.010$	$2.888\pm 0.007$	$2.890\pm 0.005$	$2.893\pm 0.003$	$2.896\pm 0.003$
Boots	$2.698\pm 0.024$	$2.822\pm 0.014$	$2.858\pm 0.010$	$2.882\pm 0.007$	$2.900\pm 0.005$	$2.905\pm 0.004$
Boots-t	$6.165\pm 0.688$	$3.832\pm 0.092$	$3.456\pm 0.040$	$3.208\pm 0.016$	$3.053\pm 0.006$	$3.007\pm 0.004$
gMult	$3.034\pm 0.194$	$2.931\pm 0.036$	$2.918\pm 0.019$	$2.912\pm 0.009$	$2.914\pm 0.005$	$2.913\pm 0.004$
gMult-t	$3.432\pm 0.045$	$3.088\pm 0.018$	$3.005\pm 0.011$	$2.951\pm 0.007$	$2.922\pm 0.005$	$2.915\pm 0.004$
rMult	$2.769\pm 0.234$	$2.758\pm 0.039$	$2.798\pm 0.018$	$2.838\pm 0.009$	$2.873\pm 0.005$	$2.888\pm 0.004$
rMult-t	$4.110\pm 0.133$	$3.350\pm 0.027$	$3.179\pm 0.015$	$3.060\pm 0.009$	$2.984\pm 0.005$	$2.960\pm 0.004$

Equations

$$\text{(E1)}\text{leavevmode}\text{leavevmode}\text{leavevmode}{\tau }_N\frac{{\hat{\eta }}_N-\eta }{\varsigma }\stackrel{N\to \infty }{\Rightarrow }\mathcal{G}(0,\mathfrak{r}),\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{(E2)}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathbb{P}\left(\underset{N\to \infty }{\lim }\Vert {\hat{\varsigma }}_N-\varsigma {\Vert }_{\infty }=0\right)=1.$$

$${\tau }_N\frac{{\hat{\eta }}_N-\eta }{{\hat{\varsigma }}_N}\stackrel{N\to \infty }{\Rightarrow }\mathcal{G}(0,\mathfrak{r}).$$

$$SCB(s,q_{\alpha,N})=\Big{[}\hat{\eta}_{N}(s)-q_{\alpha,N}\tfrac{\hat{\varsigma}_{N}(s)}{\tau_{N}},\ \hat{\eta}_{N}(s)+q_{\alpha,N}\tfrac{\hat{\varsigma}_{N}(s)}{\tau_{N}}\Big{]}$$

$$\mathbb{P}\Big{(}\forall s\in\mathcal{S}:\eta(s)\in SCB(s,q_{\alpha,N})\Big{)}=1-\alpha\,,$$

$$\mathbb{P}\Bigg{(}\max_{s\in\mathcal{S}}\tau_{N}\left|\frac{\hat{\eta}_{N}(s)-\eta(s)}{\hat{\varsigma}_{N}(s)}\right|>q_{\alpha,N}\Bigg{)}=\alpha\,.$$

$$\mathbb{P}\Big{(}\max_{s\in\mathcal{S}}\left|\mathcal{G}(0,\mathfrak{r})\right|>q_{\alpha}\Big{)}\geq\alpha\,,$$

$$Y(s)=\mu (s)+\sigma (s)Z(s),\text{ for }s\in \mathcal{S}.$$

$$\big{|}Z(s)-Z(s^{\prime})\big{|}\leq A\delta(s,s^{\prime})\text{ for all }s,s^{\prime}\in\mathcal{S}$$

$$\mathbb{E}\left[\big{(}Z(s)-Z(s^{\prime})\big{)}^{2}\right]<f(|s-s^{\prime}|)\text{ for }|s-s^{\prime}|\text{ small and }\int_{0}^{1}y^{-\frac{3}{2}}f^{\frac{1}{2}}(y)\,dy<\infty$$

$$T(s)=\frac{\sqrt{N}{G}_N(s)}{\sqrt{\frac{1}{N-1}{\sum }_{n=1}^{N-1}{G}_n(s{)}^2}},$$

$$\mathbb{E}\!\left[\big{(}D^{(d,l)}\mathcal{G}(s)-D^{(d,l)}\mathcal{G}(s^{\prime})\big{)}^{2}\right]\leq K\big{|}\log\|s-s^{\prime}\|\big{|}^{-(1+\gamma)}$$

$$\mathbb{E}\left[\chi\big{(}A(T,u)\big{)}\right]=L_{0}(\mathcal{S})\rho_{0}^{t_{N-1}}(u)+\sum_{d=1}^{D}L_{d}(\mathcal{S},\mathcal{G})\rho_{d}^{t_{N-1}}(u)\,,\leavevmode\nobreak\ D\in\mathbb{N}\,,$$

$$\frac{1}{2}\mathbb{P}\Big{(}\max_{s\in\mathcal{S}}|T(s)|>u\Big{)}\leq\mathbb{P}\Big{(}\max_{s\in\mathcal{S}}T(s)>u\Big{)}\approx\mathbb{E}\left[\chi(u)\right]\,,$$

$$\frac{1}{2}\mathbb{P}\Big{(}\max_{s\in\mathcal{S}}|T(s)|>u\Big{)}\leq\mathbb{P}\Big{(}\chi\big{(}T(s)>u\big{)}\geq 1\Big{)}\leq\mathbb{E}\left[\chi(u)\right]\,.$$

$${\widehat{\mathrm{EEC}}}_{t_{N-1}}(u)=L_0(\mathcal{S}){\rho }_0^{t_{N-1}}(u)+\sum _{d=1}^D{\hat{L}}_d(\mathcal{S},\mathcal{G}){\rho }_d^{t_{N-1}}(u)=\frac{\alpha }{2}.$$

$${\mathrm{EEC}}_{\mathcal{G}}(u)=L_0(\mathcal{S}){\rho }_0^{\mathcal{G}}(u)+\sum _{d=1}^D{L}_d(\mathcal{S},\mathcal{G}){\rho }_d^{\mathcal{G}}(u)=\frac{\alpha }{2},$$

$$\underset{k\in \{1,...,K\}}{\max }\frac{1}{\sqrt{N}}\sum _{n=1}^N{X}_{nk}$$

$$\underset{k\in \{1,...,K\}}{\max }\frac{1}{\sqrt{N}}\sum _{n=1}^N{g}_n{X}_{nk}$$

$$\underset{s\in \mathcal{S}}{\max }\sqrt{N}\left|\frac{{\hat{\eta }}_N(s)-\eta (s)}{{\hat{\varsigma }}_N(s)}\right|\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{ and }\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\underset{s\in \mathcal{S}}{\max }\frac{1}{\sqrt{N}}\left|\sum _{n=1}^N\frac{R_n^N(s)}{{\hat{\varsigma }}_N(s)}\right|$$

$$\lim_{N\rightarrow\infty}\Big{|}1-\alpha-\mathbb{P}\big{(}\leavevmode\nobreak\ \forall s\in\mathcal{S}:\leavevmode\nobreak\ \eta(s)\in SCB(s,\hat{q}_{\alpha,N})\big{)}\Big{|}\leq e^{-\big{(}\tfrac{1}{2}+\tfrac{1}{2\sigma_{c}^{2}}\big{)}\tilde{q}_{\alpha}^{2}}<e^{-\tfrac{\tilde{q}_{\alpha}^{2}}{2}}\,,$$

$$\hat{\mu}_{N}(s)=\overline{\textbf{Y}(s)}=\frac{1}{N}\sum_{n=1}^{N}Y(s)\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \hat{\sigma}_{N}^{2}(s)=\widehat{\rm var}_{N}\!\left[\!\,\textbf{Y}(s)\,\right]=\frac{1}{N-1}\sum_{n=1}^{N}\big{(}Y_{n}(s)-\bar{Y}(s)\big{)}^{2}\,,$$

$$SCB(s,{\hat{q}}_{\alpha ,N})={\hat{\mu }}_N(s)\pm {\hat{q}}_{\alpha ,N}\frac{{\hat{\sigma }}_N(s)}{\sqrt{N}}$$

$$Y(s)={\mu }_Y(s)+{\sigma }_Y(s)Z_Y(s)\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{ and }\text{leavevmode}\text{leavevmode}\text{leavevmode}X(s)={\mu }_X(s)+{\sigma }_X(s)Z_X(s),$$

$$\frac{\sqrt{N+M-2}\big{(}\overline{\textbf{Y}}-\overline{\textbf{X}}-\mu_{Y}+\mu_{X}\big{)}}{\sqrt{(1+c^{-1})\hat{\sigma}_{N}(\textbf{Y})^{2}+(1+c)\hat{\sigma}_{N}(\textbf{X})^{2}}}\xRightarrow{N,M\rightarrow\infty}\mathcal{G}=\frac{\sqrt{1+c^{-1}}\sigma_{Y}\mathcal{G}_{Y}-\sqrt{1+c}\sigma_{X}\mathcal{G}_{X}}{\sqrt{(1+c^{-1})\sigma_{Y}^{2}+(1+c)\sigma_{X}^{2}}}\,,$$

$$SCB(s,{\hat{q}}_{\alpha ,N})={\hat{\mu }}_N(s)\pm {\hat{q}}_{\alpha ,N}\frac{\sqrt{(1+c^{-1}){\sigma }_Y^2+(1+c){\sigma }_X^2}}{\sqrt{N+M}}$$

$${\Lambda }_{dl}(Z,s)={\Lambda }_{dl}(s)=\text{mbox}\mathrm{cov}\left[D^{d}Z(s),D^{l}Z(s)\right],\text{leavevmode}\text{leavevmode}\text{leavevmode}d,l=1,...,D.$$

$$L_1(\mathcal{S},Z)={\mathrm{vol}}_1(\mathcal{S},\Lambda )={\int }_S\sqrt{\text{mbox}\mathrm{var}\left[\frac{dZ}{ds}(s)\right]}ds.$$

$$\displaystyle\begin{split}L_{1}&=\frac{1}{2}{\rm length}(\partial\mathcal{S},\Lambda)=\frac{1}{2}\int_{0}^{1}\sqrt{\tfrac{d\gamma}{dt}^{T}\!\!(t)\Lambda\big{(}\gamma(t)\big{)}\tfrac{d\gamma}{dt}(t)}\,dt\\ L_{2}&={\rm vol}_{2}(\mathcal{S},\Lambda)=\int_{\mathcal{S}}\sqrt{\det\big{(}\Lambda(s)\big{)}}ds_{1}ds_{2}\,.\end{split}$$

$$R_{n}(s)=\big{(}Y_{n}(s)-\hat{\mu}(s)\big{)}/\hat{\sigma}_{N}(s)\,,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ s\in\mathcal{S}\,,n=1,...,N$$

$${\hat{L}}_1^N={\int }_0^1\sqrt{\widehat{\mathrm{var}}\left[\frac{d}{ds}\left(\frac{\text{R}(s)}{{\hat{\sigma }}_N(s)}\right)\right]}ds,$$

Full text

The theoretical work is accompanied by a simulation study comparing different methods to construct SCBs for the population mean. We show that the tGKF works well even for small sample sizes and only a Rademacher multiplier-$t$ bootstrap performs similarily well. For larger sample sizes the tGKF often outperforms the bootstrap methods and is computational faster. We apply the method to diffusion tensor imaging (DTI) fibers using a scale space approach for the difference of population means. R code is available in our Rpackage SCBfda.

1. 1 Introduction

2. 2 Simultaneous Confidence Bands

3. 2.1 SCBs for Functional Parameters

4. 2.2 Functional Signal-Plus-Noise Model

5. 3 Estimation of the Quantile

6. 3.1 Estimation of the Quantile Using the tGKF

7. 3.1.1 The Gaussian Kinematic Formula for t-processes

8. 3.1.2 The tGKF-estimator of $q_{\alpha,N}$

9. 3.2 Estimation of the Quantile Using the Bootstrap

10. 3.2.1 Non-parametric bootstrap-$t$

11. 3.2.2 Multiplier-$t$ Bootstrap

12. 4 Asymptotic Covering Rates

13. 4.1 Asymptotic SCBs for Functional Parameters

14. 4.2 Asymptotic SCBs for the Signal-Plus-Noise Model

15. 5 Estimation of the LKCs for Signal-Plus-Noise Models

16. 5.1 Definition of the LKCs for $D=1$ and $D=2$

17. 5.2 Simple Estimators Based on Residuals

18. 5.3 Estimation of LKCs in the Two Sample Problem

19. 6 Discrete Measurements, Observation Noise and Scale Spaces

20. 7 Simulations

21. 7.1 Coverage: Smooth Gaussian Case

22. 7.2 Coverage: Smooth Non-Gaussian Case

23. 7.3 Average Width and Variance of Different SCBs

24. 7.4 The Influence of Observation Noise

25. 7.5 SCBs for the Difference of Population Means of Two Independent Samples

26. 8 Application to DTI fibers

27. A Proofs

28. A.1 Proof of Proposition 1

29. A.2 Proof of Claim in Remark 3

30. A.3 Proof of Theorem 1

31. A.4 Proof of Theorem 2

32. A.5 Proofs of Theorem 3 and 4

33. A.6 Proof of Theorem 5

34. A.7 Proof of Theorem 6 and Corollary 1

35. A.8 Proof of Theorem 7

36. A.9 Proof of Corollary 2

37. A.10 Proof of Theorem 8

38. A.11 Proof of Proposition 2

1 Introduction

In the past three decades functional data analysis has received increasing interest due to the possibility of recording and storing data collected with high frequency and/or high resolution in time and space. Many different methods have been developed to study these particular data objects; for overviews of some recent developments in this fast growing field we refer the reader to the review articles Cuevas, (2014) and Wang et al., (2016).

Currently, the literature on construction of simultaneous confidence bands (SCBs) for functional parameters derived from repeated observations of functional processes is sparse. The existing methods basically split into two seperate groups. The first group is based on functional central limit theorems (fCLTs) in the Banach space of continuous functions endowed with the maximum metric and evaluation of the maximum of the limiting Gaussian process often using Monte-Carlo simulations with an estimated covariance structure of the limit process, cf. Bunea et al., (2011); Degras, (2011, 2017); Cao et al., (2012, 2014). The second group is based on the bootstrap, among others Cuevas et al., (2006); Chernozhukov et al., (2013); Chang et al., (2017), Belloni et al., (2018).

Although these methods are asymptotically achieving the correct covering probabilities, their small sample performance is often less impressive as for example discovered in Cuevas et al., (2006). Figure 1 shows the typical behaviour observed in our simulations of the covering rates of SCBs for the (smoothed) population mean using $N$ replications of a (smoothed) Gaussian process with an unknown non-isotropic covariance function. The general pattern is that SCBs using a non-parametric bootstrap-$t$ yield too wide confidence bands and thus have overcoverage, whereas a functional Gaussian multiplier bootstrap inspired by Chernozhukov et al., (2013) and Degras, (2011) asymptotic SCBs lead to undercoverage. We used here SCBs for the smoothed population mean, since the R-package SCBmeanfd by Degras does require smoothing.

Contrary to most of the current methods, which assume that each functional process is observed only at discrete points of their domain, we start from the viewpoint that the whole functional processes are observed and sufficiently smooth, since often the first step in data analysis is anyhow smoothing the raw data. This implies that we construct SCBs only for the smoothed population parameter and thereby circumventing the bias problem altogether to isolate the covering aspect. In principle, however, it is possible to apply our proposed method also to discretely observed data and inference on the true (non-smoothed) mean as for example in settings described in Zhang et al., (2007); Degras, (2011); Zhang et al., (2016). But for the sake of simplicity we restrict ourselves to only robustify against the choice of a smoothing parameter by introducing SCBs for scale space surfaces, which were proposed in Chaudhuri and Marron, (1999).

The main contribution of this paper is proposing the construction of SCBs for general functional parameters based on the Gaussian kinematic formula for $t$-processes (tGKF) and studying its theoretical covering properties. Moreover, as an important case study we explain and prove how all the required assumptions can be satisfied for construction of SCBs for the population mean in functional signal-plus-noise models and scale spaces. A first glance on the improvement, especially for small samples, using the tGKF is shown in Figure 1.

A second contribution is that we compare to a studentized version of the multiplier bootstrap based on residuals, which we call multiplier-$t$ bootstrap (Mult-$t$), which also improves the small sample properties of the bootstrap SCBs considerably, if the multiplier is chosen to be Rademacher random variables. Since the scope of the paper is to analyse the theoretical properties of the use of the tGKF in construction of SCBs, we do not provide any rigorous mathematical theory about the Mult-$t$ bootstrap. Theoretical analysis of this bootstrap method, in particular the influence of the choice of multiplier, is an interesting opportunity for future research. For recent work on the performance of different multipliers in similar scenarios, see e.g., Deng and Zhang, (2017) and Davidson and Flachaire, (2008).

SCBs usually require estimation of the quantiles of the maximum of a process. Often this process is the maximum of a process derived from a functional central limit theorem (fCLT); thus the maximum of a Gaussian process. Our proposed method tries to resolve the poor covering rates for small sample sizes by imposing two small but significant changes. Firstly, we use a pointwise $t$-statistic, since for small sample sizes and unknown variance this would be the best choice to construct confidence bands for Gaussian random variables. Secondly, we employ a formula which is known to approximate the tail distribution of the maximum of a pointwise $t$-distributed process very well. This formula is known as the Gaussian kinematic Formula for $t$-processes (tGKF).

In a nutshell the tGKF as proven in even more generality in Taylor, (2006) expresses the expected Euler characteristic (EEC) of the excursion set of a random process $F(Z_{1},...,Z_{N})$, $F\in\mathcal{C}^{2}(\mathbb{R}^{N},\mathbb{R})$, derived from unit variance Gaussian processes $Z_{1},...,Z_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}Z$ indexed by a nice compact subset of $\mathbb{R}^{D}$ in terms of a finite sum of $N$ known functions with corresponding $N$ coefficients. These coefficients are called Lipschitz-Killing curvatures (LKCs) and depend solely on the random process $Z$. Earlier versions of this formula derived in Adler, (1981) for Gaussian processes were used by Worsley et al., (1996),Worsley et al., (2004) in neuroscience for multiple testing corrections. Takemura and Kuriki, (2002) have shown that the GKF (for Gaussian processes) is closely related to the Volume of Tubes formula dating all the way back to Working and Hotelling, (1929) and which has been applied for SCBs in nonlinear regression analysis, e.g., Johansen and Johnstone, (1990); Krivobokova et al., (2010); Lu and Kuriki, (2017).

In this sense the version of the tGKF of Taylor, (2006) can be interpreted as a generalization of the volume of tube formula for repeated observations of functional data. The most important and valuable aspect of the tGKF is that it is a non-asymptotic formula in the sample size $N$. Hence it is suitable for applications with small sample sizes. Moreover, it has been proven in Taylor et al., (2005) that the EEC is a good approximation of the quantile of the maximum value of Gaussian- and $\chi^{2}$-processs over their domain. Although it is not yet proven, simulations suggest that this property seems to be also true for more complicated Gaussian related processes like $t$-processes. Therefore, it can be used to control the family wise error rate (FWER) as in Taylor and Worsley, (2007); Taylor et al., (2009). To the best of our knowledge, the full potential of the tGKF for functional data seems to not have been explored yet in the statistical literature and we will try to fill this gap and tie together some loose strings for the challenge of constructing SCBs.

Our main theoretical contributions are the following. In Theorem 2 we show based on the main result in Taylor et al., (2005) that asymptotically the error in the covering rate of SCBs for a function-valued population parameter based on the GKF for $t$-processes can be bounded and is small, if the targeted covering probability of the SCB is sufficiently high. This requires no Gaussianity of the observed processes. It only requires that the estimator of the targeted function-valued parameter fulfills a fCLT in the Banach space of continuous functions with a sufficiently regular Gaussian limit process. Moreover, it requires consistent estimators for the LKCs. Using this general result we derive SCBs for the population mean and the difference of population means for functional signal-plus-noise models, where we allow the error processes to be non-Gaussian. Especially we derive for such models defined over sufficiently regular domains $\mathcal{S}\subset\mathbb{R}^{D}$, $D=1,2$, consistent estimators for the LKCs. These estimators are closely related to the discrete estimators in Taylor and Worsley, (2007), but we can even prove CLTs for our estimators. In order to deal with observation noise we give in Theorem 8 sufficient conditions to have weak convergence of a scale space process to a Gaussian limit extending the results from Chaudhuri and Marron, (2000) from regression analysis to repeated observations of functional data. Additionally, we prove that also the LKCs of this limit process can be consistently estimated and therefore Theorem 2 can be used to bound the error in the covering rate for SCBs of the population mean of a scale space process.

These theoretical findings are accompanied by a simulation study using our Rpackage SCBfda, which can be found on https://github.com/ftelschow/SCBfda and a data application to diffusion tensor imaging (DTI) fibers using a scale space approach for the difference of population means. In the simulation study we compare the performance of the tGKF approach to SCBs for different error processes mainly with bootstrap approaches and conclude that the tGKF approach does not only often give better coverings for small sample sizes, but also outperforms bootstrap approaches computationally. Moreover, the average width of the tGKF confidence bands is lower for large sample sizes.

Organisation of the Article

Our article is organized as follows: In Section 2.1 we describe the general idea of construction of SCBs for functional parameters. Section 2.2 defines the functional signal-plus-noise model and general properties on the error processes, which are necessary to prove most of our theorems. The next section explains how SCBs can be constructed in praxis and especially explains how the tGKF can be used for this purpose. Section 4 finally proves asymptotic properties of the SCBs constructed using the tGKF for general functional parameters. The latter require consistent estimation of the LKCs, which are discussed for the functional signal-plus-noise model in Section 5 together with the new CLTs for the LKCs. Robustification using SCBs for Scale Space models can be found in 6. In Section 7 we compare our proposed method in different simulations to competing methods to construct SCBs in the case of the functional signal-plus-noise model for various settings, which is followed in Section 8 by a data example.

2 Simultaneous Confidence Bands

2.1 SCBs for Functional Parameters

We describe now a general scheme for construction of simultaneous confidence bands (SCBs) for a functional parameter $s\mapsto\eta(s)$, $s\in\mathcal{S}$, where $\mathcal{S}$ is a compact metric space. Note that all functions of $s$ hereafter will be assumed to belong to $\mathcal{C}(\mathcal{S})$ the space of continuous functions from $\mathcal{S}$ to $\mathbb{R}$.

Assume we have estimators $s\mapsto\hat{\eta}_{N}(s)$ of $\eta$ and $s\mapsto\hat{\varsigma}_{N}(s)$ estimating a different functional parameter $s\mapsto\varsigma(s)$ fulfilling

[TABLE]

Here and whenever convergence of random continuous functions is considered, “$\Rightarrow$” denotes weak convergence in $\mathcal{C}(\mathcal{S})$ endowed with the maximums norm $\|\cdot\|_{\infty}$. In (E1), $\mathcal{G}(0,\mathfrak{r})$ is a mean zero Gaussian process with covariance function $\mathfrak{r}$ satisfying $\mathfrak{r}(s,s)=1$ for all $s\in\mathcal{S}$, $\{\tau_{N}\}_{N\in\mathbb{N}}$ is a scaling sequence of positive numbers. Assumptions (E1) and (E2) together with the functional version of Slutzky’s Lemma imply

[TABLE]

Thus, it is easy to check that the collection of intervals

[TABLE]

form $(1-\alpha)$-simultaneous confidence bands of $\eta$, i.e.

[TABLE]

provided that

[TABLE]

Unfortunately, the quantiles $q_{\alpha,N}$ are in general unknown and need to be estimated. Therefore, the main contribution of this article is the study and comparison of different estimators for $q_{\alpha,N}$.

In principle, there are two general approaches. Limit approximations try to estimate $q_{\alpha,N}$ by approximations of the quantiles of the asymptotic process, i.e.

[TABLE]

as done in Degras, (2011, 2017) for the special case of a signal-plus-noise model and the local linear estimator of the signal. Here usually the covariance function $\mathfrak{r}$ is estimated and then many samples of the Gaussian process $\mathcal{G}(0,\hat{\mathfrak{r}})$ are simulated in order to approximate $q_{\alpha}$.

Better performance for finite sample sizes can be achieved by approximating $q_{\alpha,N}$ directly by bootstrap approaches such as a fully non-parametric bootstrap proposed in Degras, (2011), which is inspired by the bootstrap-$t$ confidence intervals (e.g., DiCiccio and Efron, (1996)),.

Moreover, the tGKF approach, which we will introduce in Section 3.1, uses a “parametric” estimator of $q_{\alpha,N}$, approximating the l.h.s of (1) by a $t$-process over $\mathcal{S}$.

2.2 Functional Signal-Plus-Noise Model

Our leading example will be SCBs for the population mean curve in a functional signal-plus-noise model which we will introduce now. Let us assume that $\mathcal{S}\subset\mathbb{R}^{D}$, $D\in\mathbb{N}$, is a compact set with piecewise $\mathcal{C}^{2}$-boundary $\partial\mathcal{S}$. The functional signal-plus-noise model is given by

[TABLE]

Here we assume that $\mu,\sigma$ are continously differentiable functions on $\mathcal{S}$ and $Z$ is a stochastic process with zero mean and covariance function $\mbox{\rm cov}\left[Z(s),Z(s^{\prime})\right]=\mathfrak{c}(s,s^{\prime})$ for $s,s^{\prime}\in\mathcal{S}$ satisfying $\mathfrak{c}(s,t)=1$ if and only if $s=t$. Moreover, we introduce the following useful properties of stochastic processes.

Definition 1**.**

We say a stochastic process $Z$ with domain $\mathcal{S}$ is $(\mathcal{L}^{p},\delta)-$Lipschitz, if there is a (semi)-metric $\delta$ on $\mathcal{S}$ continuous w.r.t. the standard metric on $\mathbb{R}^{D}$ and a random variable $A$ satisfying \mathbb{E}\big{[}|A|^{p}\big{]}<\infty such that

[TABLE]

and $\int_{0}^{1}H^{1/2}(\mathcal{S},\delta,\epsilon)d\epsilon<\infty$, where $H(\mathcal{S},\delta,\epsilon)$ denotes the metric entropy function of the (semi)-metric space $(\mathcal{S},\delta)$, e.g., Adler and Taylor, (2007, Def. 1.3.1.).

Remark 1**.**

Any $(\mathcal{L}^{2},\delta)-$Lipschitz process $Z$ has necessarily almost surely continuous sample paths. Moreover, this property is the main ingredient in the version of a CLT in $\mathcal{C}(\mathcal{S})$ proven in Jain and Marcus, (1975) or Ledoux and Talagrand, (2013, Section 10.1). However, there are different results on CLTs in $\mathcal{C}(\mathcal{S})$ with different assumptions on the process, which in principle could replace this condition. For example for $D=1$ we could also use the condition

[TABLE]

where $f:[0,1]\rightarrow\mathbb{R}_{\geq 0}$ is non-decreasing near 0 and satisfies $f(0)=0$. This is due to Hahn, (1977). However, $(\mathcal{L}^{2},\delta)$-Lipschitz seems to be the most tractable assumption for our purposes.

Definition 2**.**

We say a process $Z$ has finite $p$-th $\mathcal{C}(\mathcal{S})$-moment, if \mathbb{E}\big{[}\|Z(s)\|^{p}_{\infty}\big{]}<\infty.

Proposition 1**.**

Any $(\mathcal{L}^{p},\delta)-$Lipschitz process over a compact set $\mathcal{S}$ has finite $p$-th $\mathcal{C}(\mathcal{S})$-moment.

Remark 2**.**

Since any continuous Gaussian process satisfies the finite $p$-th $\mathcal{C}(\mathcal{S})$-moment condition, cf. Landau and Shepp, (1970), it is possible to proof a reverse of Proposition 1 for continuously differentiable Gaussian processes. Moreover, finite $p$-th $\mathcal{C}(\mathcal{S})$-moment conditions are often assumed, when uniform consistency of estimates of the covariance function are required, e.g., Hall et al., (2006) and Li et al., (2010).

3 Estimation of the Quantile

This section describes different estimators for the quantiles $q_{\alpha,N}$ defined by equation (3). Especially, we propose using the Gaussian kinematic formula for $t$-processes (tGKF) as proven in Taylor, (2006). Moreover, we describe different bootstrap estimators.

3.1 Estimation of the Quantile Using the tGKF

3.1.1 The Gaussian Kinematic Formula for t-processes

A $t_{N-1}$-process $T$ over an index set $\mathcal{S}$ is a stochastic process such that $T(s)$ is $t_{N-1}$-distributed for all $s\in\mathcal{S}$. In order to match the setting in Taylor, (2006) we assume that the process is given by

[TABLE]

where $G_{1},...,G_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}\mathcal{G}$ are mean zero, variance one Gaussian processes. Let us define $D^{I}\!\mathcal{G}=\tfrac{\partial^{|I|}\mathcal{G}}{\partial s_{I_{1}}...\partial s_{I_{K}}}$, where $K=|I|$ denotes the number of elements in the multi-index $I$. Then we require the following assumptions:

1. (G1)

$\mathcal{G}$ has almost surely $\mathcal{C}^{2}-$sample paths. 2. (G2)

The joint distribution of \Big{(}D^{d}\mathcal{G}(s),D^{(d,l)}\mathcal{G}(s)\Big{)} is nondegenerate for all $s\in\mathcal{S}$ and $d,l=1,...,D$. 3. (G3)

There is an $\epsilon>0$ such that

[TABLE]

for all $d,l=1,...,D$ and for all $|s-t|<\epsilon$. Here $K>0$ and $\gamma>0$ are finite constants.

Remark 3**.**

Assumption (G3) is satisfied for any process $\mathcal{G}$ having almost surely $\mathcal{C}^{3}-$sample paths and all third derivatives have finite second $\mathcal{C}(\mathcal{S})$-moment, see Definiton 2. In particular, this holds true for any Gaussian process with $\mathcal{C}^{3}-$sample paths. For completeness the argument is carried out in more detail in the appendix.

Under these assumptions the tGKF is an exact, analytical formula of the expectation of the Euler characteristic $\chi$ of the excursion sets $A(T,u)=\left\{s\in\mathcal{S}\leavevmode\nobreak\ |\leavevmode\nobreak\ T(s)>u\right\}$ of $T$. This formula as proven in Taylor, (2006) or Adler and Taylor, (2007) is

[TABLE]

where $\rho_{d}^{t_{N-1}}$, $d=1,...,D$, is the $d$-th Euler characteristic density of a $t_{N-1}$-process, which can be found for example in Taylor and Worsley, (2007). Moreover, $L_{d}$ denotes the $d$-th Lipschitz killing curvature, which only depends on $\mathcal{G}$ and the parameter space $\mathcal{S}$. Note that $L_{0}(\mathcal{S})=\chi(\mathcal{S})$.

Equation (8) is useful, since by the expected Euler characteristic heuristic (EECH) (see Taylor et al., (2005) ), we expect

[TABLE]

to be a good approximation for large thresholds $u$. In the case that $\mathcal{S}\subset\mathbb{R}$ this approximation actually is always from above. This is due to the fact that the Euler characteristic of one-dimensional sets is always non-negative and hence using the Markov-inequality we obtain

[TABLE]

The same argument is heuristically valid for high enough thresholds in any dimension, since the excursion set will with high probabilty consist mostly of simply-connected sets. Thus, the EC of the excursion set will be positive with high probability. Therefore, once again we expect that the EEC is an approximation from above for the excursion probability.

Notably, the tGKF Equation (8) had a predecessor for the EEC of a Gaussian processes $\mathcal{G}$. Remarkably, the only difference between these formulas is that the Euler characteristic densities are different, see Adler and Taylor, (2007, p.315, (12.4.2)). Moreover, the approximation error when using the EEC instead of the excursion probabilities of the Gaussian process $\mathcal{G}$ has been analytically quantified in Taylor et al., (2005).

3.1.2 The tGKF-estimator of $q_{\alpha,N}$

Assume that we have consistent estimators $\hat{L}_{d}(\mathcal{S},\mathcal{G})$ for the LKCs $L_{d}(\mathcal{S},\mathcal{G})$ for all $d=1,...,D$ (we provide such estimators for $D=1,2$ in Section 5). Then a natural estimator $\hat{q}_{\alpha,N}$ of the quantile $q_{\alpha,N}$ is the largest solution $u$ of

[TABLE]

The following result will be used in the next section for the proof of the accuracy of the SCBs derived using the tGKF estimator $q_{\alpha,N}$.

Theorem 1**.**

Assume that $\hat{L}^{N}_{d}(\mathcal{S},\mathcal{G})$ is a consistent estimator of $L_{d}(\mathcal{S},\mathcal{G})$. Then $\hat{q}_{\alpha,N}$ given by (11) converges almost surely for $N$ tending to infinity to the largest solution $\tilde{q}_{\alpha}$ of

[TABLE]

where $\rho_{D}^{\mathcal{G}}$ are the Euler characteristic densities of a Gaussian process, which can be found in Adler and Taylor, (2007, p.315, (12.4.2)).

3.2 Estimation of the Quantile Using the Bootstrap

As alternatives to the approach using the tGKF we discuss a non-parametric bootstrap-$t$ and a multiplier bootstrap estimator of the quantile $q_{\alpha,N}$ defined in Equation (3).

3.2.1 Non-parametric bootstrap-$t$

Based on Degras, (2011) we review a bootstrap estimator of the quantile $q_{\alpha,N}$. Assume that the estimators $s\mapsto\hat{\eta}_{N}(s)$ and $s\mapsto\hat{\varsigma}_{N}(s)$ are obtained from a sample $Y_{1},...,Y_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}Y$ of random functions, then the non-parametric bootstrap-$t$ estimator of $q_{\alpha,N}$ is obtained as follows

Resample from $Y_{1},...,Y_{N}$ with replacement to produce a bootstrap sample $Y_{1}^{*},...,Y_{N}^{*}$. 2. 2.

Compute $\hat{\eta}^{*}_{N}$ and $\hat{\varsigma}_{N}^{*}$ using the sample $Y_{1}^{*},...,Y_{N}^{*}$. 3. 3.

Compute T^{*}=\max_{s\in\mathcal{S}}\tau_{N}\big{|}\hat{\eta}^{*}_{N}(s)-\hat{\eta}_{N}(s)\big{|}/\hat{\varsigma}_{N}^{*}(s). 4. 4.

Repeat steps 1 to 3 many times to approximate the conditional law \mathcal{L}^{*}=\mathcal{L}\big{(}T^{*}\,|\leavevmode\nobreak\ Y_{1},...,Y_{N}\big{)} and take the $(1-\alpha)\cdot 100\%$ quantile of $\mathcal{L}^{*}$ to estimate $q_{\alpha,N}$.

Remark 4**.**

Note that the variance in the denominator is also bootstrapped, which corresponds to the standard bootstrap-$t$ approach for confidence intervals, cf. DiCiccio and Efron, (1996). This is done in order to mimic the l.h.s. in (1), and improves the small sample coverage. Results of this bootstrap variant will be discussed in our simulations.

We expect that this estimator works well for large enough sample sizes. Although Degras, (2011) introduced it especially for small sample sizes, there is not much hope that it will perform well in this case, since it is well known that confidence bands for a finite dimensional parameter based on the bootstrap-$t$ have highly variable end points, cf., Good, (2005, Section 3.3.3). Evidence that this remains the case in the functional world will be given in our simulations in Section 7.

3.2.2 Multiplier-$t$ Bootstrap

The second bootstrap method, which we introduce, builds on residuals and a version of the multiplier (or Wild) bootstrap (e.g., Mammen, (1993)) designed for the maximum of sums of $N$ independent random variables in high dimensions as discussed in detail by Chernozhukov et al., (2013). Briefly, the approach in Chernozhukov et al., (2013) is as follows. Let $X_{1},...,X_{N}$ be independent random vectors in $\mathbb{R}^{K}$, $N,K\in\mathbb{N}$ with \mathbb{E}\big{[}X_{n}\big{]}=0 and finite covariance \mathbb{E}\big{[}X_{n}X_{n}^{T}\big{]} for all $n\in\{1,...,N\}$. Define X_{n}^{T}=\big{(}X_{n1},...,X_{nK}\big{)} and assume there are $c,C\in\mathbb{R}_{>0}$ such that c<\mathbb{E}\big{[}X_{nk}^{2}\big{]}<C for all $n\in\{1,...,N\}$ and all $k\in\{1,...,K\}$. Under these assumptions it is shown in Chernozhukov et al., (2013, Theorem 3.1) that the quantiles of the distribution of

[TABLE]

can be asymptotically consistently estimated by the quantiles of the multiplier bootstrap i.e., by the distribution of

[TABLE]

with multipliers $g_{1},...,g_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}\mathcal{N}(0,1)$ given the data $X_{1},...,X_{N}$, even if $K\gg N$.

We adapt this approach to the estimation of $q_{\alpha,N}$ for functional parameters. Therefore we again assume that the estimators $s\mapsto\hat{\eta}_{N}(s)$ and $s\mapsto\hat{\varsigma}_{N}(s)$ are obtained from an i.i.d. sample $Y_{1},...,Y_{N}$ of random functions. However, here we also need to assume that we have residuals $R_{n}^{N}$ for $n=1,...,N$ satisfying $R_{n}^{N}/\varsigma\xRightarrow{N\rightarrow\infty}\mathcal{G}(0,\mathfrak{r})$ and mutual independence for $N$ tending to infinity. For example for the signal-plus-noise model, where $\hat{\eta}_{N}=\hat{\mu}_{N}$ and $\hat{\varsigma}_{N}=\hat{\sigma}_{N}$ are the pointwise sample mean and the pointwise sample standard deviation, the residuals $R_{n}^{N}=\sqrt{\tfrac{N}{N-1}}(Y_{n}-\hat{\mu}_{N})$ do satisfy these conditions, if the error process $Z$ is $(\mathcal{L}^{2},\delta)$-Lipschitz and has finite second $\mathcal{C}(\mathcal{S})$-moment. Thus, using these assumptions we obtain that

[TABLE]

have approximately the same distribution. Algorithmitically, our proposed multiplier bootstrap is

Compute residuals $R_{1}^{N},...,R_{N}^{N}$ and multipliers $g_{1},...,g_{N}\overset{\text{i.i.d.}}{\sim}g$ with $\mathbb{E}[g]=0$ and $\mbox{\rm var}[g]=1$ 2. 2.

Estimate $\hat{\varsigma}_{N}^{*}(s)$ from $g_{1}Y_{1}(s),...,g_{N}Y_{N}(s)$. 3. 3.

Compute $T^{*}(s)=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}g_{n}\frac{R_{n}^{N}\!(s)}{\hat{\varsigma}_{N}^{*}(s)}$. 4. 4.

Repeat steps 1 to 3 many times to approximate the conditional law \mathcal{L}^{*}=\mathcal{L}\big{(}T^{*}\,|\leavevmode\nobreak\ Y_{1},...,Y_{N}\big{)} and take the $(1-\alpha)\cdot 100\%$ quantile of $\mathcal{L}^{*}$ to estimate $q_{\alpha,N}$.

In our simulations we use Gaussian multipliers as proposed in Chernozhukov et al., (2013), but find that Rademacher multipliers as used in regression models with heteroskedastic noise, e.g. Davidson and Flachaire, (2008), perform much better for small sample sizes.

Remark 5**.**

Again the choice to compute bootstrapped versions of $\hat{\varsigma}_{N}$ mimics the true distribution of the maximum of the random process on the l.h.s of (1) for finite $N$ better than just applying the multipliers to $R_{n}^{N}/\hat{\varsigma}_{N}$.

4 Asymptotic Covering Rates

4.1 Asymptotic SCBs for Functional Parameters

This section discusses the accuracy of the SCBs derived using the tGKF. Since the expected Euler characteristic of the excursion sets are only approximations of the excursion probabilities, there is no hope to prove that the covering of these confidence bands is exact. Especially, if $\alpha$ is large the approximation fails badly and will usually lead to confidence bands that are too wide. However, for values of $\alpha<0.1$ typically used for confidence bands, the EEC approximation works astonishingly well. Theoretically, this has been made precise for Gaussian processes in Theorem 4.3 of Taylor et al., (2005), which is the main ingredient in the proof of the next result. Additionally, it relies on the fCLT (1) and the consistency of $q_{\alpha,N}$ for $q_{\alpha}$ proved in Theorem 1.

Theorem 2**.**

Assume (E1-2) and assume that the limiting Gaussian process $\mathcal{G}(0,\mathfrak{r})$ satisfies (G1-3). Moreover, let $\hat{L}_{d}^{N}$ be a sequence of consistent estimators of $L_{d}$ for $d=1,...,D$. Then there exists an $\alpha^{\prime}\in(0,1)$ such that for all $\alpha\leq\alpha^{\prime}$ we have that the SCBs defined in Equation (2) fullfill

[TABLE]

where $\sigma_{c}^{2}$ is the critical variance of an associated process of $\mathcal{G}(0,\mathfrak{r})$, $\hat{q}_{\alpha,N}$ is the quantile estimated using equation (11) and $\tilde{q}_{\alpha}$ is defined in Theorem 1.

Typically, in our simulations we have that, for $\alpha=0.05$, the quantile $\tilde{q}_{\alpha}$ is about $3$ leading to an upper bound of $\approx 0.011$, if we use the weaker bound without the critical variance.

4.2 Asymptotic SCBs for the Signal-Plus-Noise Model

As an immediate application of Theorem 2 we derive now SCBs for the population mean and the difference of population means in one and two sample scenarios of the signal-plus-noise model introduced in Section 2.2. The derivation of consistent estimators for the LKCs will be postponed to the next section.

Theorem 3** (Asymptotic SCBs for Signal-Plus-Noise Model).**

Let $Y_{1},...,Y_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}Y$ be a sample of model (5) and assume $Z$ is an $(\mathcal{L}^{2},\delta)-$Lipschitz process. Define \textbf{Y}(s)=\big{(}Y_{1}(s),...,Y_{N}(s)\big{)}.

1. (i)

Then the estimators

[TABLE]

fullfill the conditions (E1-2) with $\tau_{N}=\sqrt{N}$, $\eta=\mu$, $\varsigma=\sigma$ and $\mathfrak{r}=\mathfrak{c}$. 2. (ii)

*If there are consistent estimators of the LKCs $L_{d}$ and $Z$ has $\mathcal{C}^{3}$-sample paths and all partial derivatives up to order $3$ of $Z$ are *$(\mathcal{L}^{2},\delta)-$Lipschitz processes with finite $\mathcal{C}(\mathcal{S})$-variances and $\mathcal{G}(0,\mathfrak{c})$ fullfills the non-degeneracy condition (G2), then the accuracy result of Theorem 2 holds true for the SCBs

[TABLE]

with $\hat{q}_{\alpha,N}$ of Theorem 1.

Remark 6**.**

A simple condition on $Z$ to ensure that $\mathcal{G}(0,\mathfrak{c})$ fulfills the non-degeneracy condition (G2) is that for all $d,l\in\{1,...,D\}$ we have that $\mbox{\rm cov}\!\left[\left(D^{l}Z(s),D^{(d,l)}Z(s)\right)\right]$ has full rank for all $s\in\mathcal{S}$. A proof is provided in Lemma 5 in the appendix.

Theorem 4** (Asymptotic SCBs for Difference of Means of Two Signal-Plus-Noise Models).**

Let $Y_{1},...,Y_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}Y$ and $X_{1},...,X_{M}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}X$ be independent samples, where

[TABLE]

with $Z_{Y},Z_{X}$ both $(\mathcal{L}^{2},\delta)-$Lipschitz processes and assume that $c=\lim_{N,M\rightarrow\infty}N/M$. Then

1. (i)

Condition (E1) is satisfied, i.e.

[TABLE]

where $G_{Y},G_{X}$ are Gaussian processes with the same covariance structures as $X$ and $Y$ and the denominator converges uniformly almost surely, i.e. condition (E2) is satisfied. 2. (ii)

*If there are consistent estimators of the LKCs $L_{d}$ of $\leavevmode\nobreak\ \mathcal{G}$ and $Z_{X},Z_{Y}$ have $\mathcal{C}^{3}$-sample paths, fullfill the non-degeneracy condition (G2) and all their partial derivatives are *$(\mathcal{L}^{2},\delta)-$Lipschitz processes with finite $\mathcal{C}(\mathcal{S})$-variances, then the accuracy result of Theorem 2 holds true for the SCBs

[TABLE]

with $\hat{q}_{\alpha,N}$ of Theorem 1.

5 Estimation of the LKCs for Signal-Plus-Noise Models

We turn now to consistent estimators of the LKCs. However, we treat only the case where $Z$ is a process defined over a domain $\mathcal{S}$, where $\mathcal{S}$ is assumed to be either a compact collection of intervals in $\mathbb{R}$ or a compact two-dimensional domain of $\mathbb{R}^{2}$ with a piecewise $\mathcal{C}^{2}$-boundary $\partial\mathcal{S}$. We restrict ourselves to this setting, since there are simple, explicit formulas for the LKCs.

5.1 Definition of the LKCs for $D=1$ and $D=2$

The basic idea behind the Lipschitz-Killing curvatures $L_{d}(\mathcal{S},Z)$ is that they are the intrinsic volumes of $\mathcal{S}$ with respect to the pseudo-metric $\tau(s,s^{\prime})=\sqrt{\mbox{\rm var}\left[Z(s)-Z(s^{\prime})\right]}\,,$ which is induced by the pseudo Riemannian metric given in standard coordinates of $\mathbb{R}^{D}$, cf. Taylor and Adler, (2003),

[TABLE]

For the cases $D=1,2$ the general expressions from Adler and Taylor, (2007) can be nicely simplified. For $D=1$, i.e. $\mathcal{S}\subset\mathbb{R}$, the only unknown is $L_{1}$, which is given by

[TABLE]

In the case of $\mathcal{S}\subset\mathbb{R}^{2}$ with piecewise $\mathcal{C}^{2}$-boundary $\partial\mathcal{S}$ parametrized by the piecewise $\mathcal{C}^{2}$-function $\gamma:[0,1]\rightarrow\mathbb{R}^{2}$ the LKCs are given by

[TABLE]

Note that by the invariance of the length of a curve the particular choice of the parametrization $\gamma$ of $\mathcal{S}$ does not change the value of $L_{1}$.

5.2 Simple Estimators Based on Residuals

In order to estimate the unknown LKCs $L_{1}$ and $L_{2}$ from a sample $Y_{1},...,Y_{N}\sim Y$ of model (5), we assume that the estimators $\hat{\mu}_{N}$ of $\mu$ and $\hat{\sigma}_{N}$ of $\sigma$ satisfy (E1-2) with $\mathfrak{r}=\mathfrak{c}$ and

1. (L)

$\lim_{N\rightarrow\infty}\|D^{d}\!\left(\sigma/\hat{\sigma}_{N}\right)\|_{\infty}=0$ for all $d=1,...,D$ almost surely.

These conditions imply that the normed residuals

[TABLE]

and its gradient converge uniformly almost surely to the random processes $Z_{n}$ and $\nabla Z_{n}$. Let us denote the vector of residuals by \textbf{R}(s)=\big{(}R_{1}(s),...,R_{N}(s)\big{)}. In view of equation (16) it is hence natural to estimate the LKC $L_{1}$ for $D=1$ by

[TABLE]

and using equations (17) the LKCs for $D=2$ by

[TABLE]

where $\hat{\Lambda}_{N}(s)=\widehat{\mbox{\rm var}}_{N}\left[\nabla\textbf{R}(s)\right]$ is the empirical covariance matrix of $\nabla\textbf{R}(s)$. Thus, we need to study mainly the properties of this estimator to ensure properties of the estimated LKCs.

Theorem 5**.**

Assume $Z$ and $D^{d}Z$ for $d=1,...,D$ are $(\mathcal{L}^{2},\delta)-$Lipschitz processes and (E2) and (L) hold true. Then, we obtain

[TABLE]

If $Z$ and $D^{d}Z$ for $d=1,...,D$ are even $(\mathcal{L}^{4},\delta)-$Lipschitz processes, we obtain

[TABLE]

with $\iota:{\rm Sym(2)}\rightarrow\mathbb{R}^{3}$ mapping

[TABLE]

and the matrix valued covariance function $\mathfrak{t}:\mathcal{S}\times\mathcal{S}\rightarrow{\rm Sym(3)}$ is given componentwise by

[TABLE]

for all $s,s^{\prime}\in\mathcal{S}$ and $d,l=1,...,D$.

From this theorem we can draw two interesting results for $D<3$. The first one is a consistency results of the LKCs. In a sense this is just reproving the result for the discrete estimator given in Taylor and Worsley, (2007) without the need to refer to convergences of meshes and making underlying assumptions more visible. Especially, it becomes obvious that the Gaussianity assumption on $Z$ is not required, since we only have to estimate the covariance matrix of the gradient consistently.

Theorem 6** (Consistency of LKCs).**

Under the setting of Theorem 5 it follows for $d=1,2$ that

[TABLE]

An immediate corollary can be drawn for the generic pointwise estimators in the functional signal-plus-noise Model.

Corollary 1**.**

Let $\hat{\sigma}^{2}_{N}(s)=\widehat{\mbox{\rm var}}_{N}\left[\textbf{Y}(s)\right]$ and assume that $Z$ and $D^{d}\!Z$ for $d=1,2$ are $(\mathcal{L}^{2},\delta)-$Lipschitz processes. Then for $d=1,2$ we obtain

[TABLE]

Remark 7**.**

If the process $Z$ is Gaussian then we can even infer that the estimators (19) and (20) are unbiased, if $\hat{\sigma}^{2}_{N}(s)=\widehat{\mbox{\rm var}}_{N}\left[\textbf{Y}(s)\right]$, see Taylor and Worsley, (2007).

The second conclusion we can draw from this approach is that we can derive a CLT for the estimator of the LKCs.

Theorem 7** (CLT for LKCs).**

Assume $Z$ and $D^{d}\!Z$ for $d=1,2$ are even $(\mathcal{L}^{4},\delta)-$Lipschitz processes and (E2) and (L) hold true. Then, if $\mathcal{S}\subset\mathbb{R}$,

[TABLE]

and, if $\mathcal{S}\subset\mathbb{R}^{2}$, we obtain

[TABLE]

where $\mathcal{G}(s)$ is a Gaussian process with mean zero and covariance function $\mathfrak{t}$ as in Theorem 5 and $\gamma$ is a parametrization of the boundary $\partial\mathcal{S}$.

Corollary 2**.**

Assume additionally to the Assumptions of Theorem 7 that $Z$ is Gaussian with covariance function $\mathfrak{c}$ and $\mathcal{S}\subset\mathbb{R}$. Then, we have the simplified representation

[TABLE]

where

[TABLE]

5.3 Estimation of LKCs in the Two Sample Problem

In the previous section we dealt with the case of estimating the LKCs in the case of one sample. This idea can be extended to the two sample case leading to a consistent estimator of the LKCs of the asymptotic process given in Theorem 4. Using the same notations as before we have that the difference in mean satisfies the following fCLT

[TABLE]

Therefore, note that independence of the samples implies that the covariance matrix defined in equation (15) splits into a sum of covariance matrices depeding on $\mathcal{G}_{X}$ and $\mathcal{G}_{Y}$, i.e.,

[TABLE]

The latter summands, however, under the assumptions of Theorem 5 can be separately consistently estimated using the previously discussed method and the residuals

[TABLE]

The sum of these estimators is a consistent estimator of $\Lambda(\mathcal{G})$ and thus the estimator described in the previous section based on the estimate of $\Lambda(\mathcal{G})$ using the above residuals gives a consistent estimator of the LKCs of $\mathcal{G}$.

6 Discrete Measurements, Observation Noise and Scale Spaces

We discuss now the realistic scenario, where the data is observed on a measurement grid with possible additive measurement noise. Again for simplicity we restrict ourselves to the signal-plus-noise model, which is

[TABLE]

where $\mathbf{S}=(s_{1},...,s_{p})\in\mathcal{S}\subset\mathbb{R}^{D\times P}$. Here we assume that $\varepsilon$ is stochastic process on $\mathcal{S}$ with finite second $\mathcal{C}(\mathcal{S})$-moment representing the observation noise and covariance function $\mathfrak{e}$. We also assume that $\varepsilon$, $Z$ and $\mathbf{S}$ are mutually independent. We are interested in doing inference on the mean $\mu$, if we observe a sample $(\mathbf{S}_{1},Y_{1}),...,(\mathbf{S}_{N},Y_{N}){\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}(\mathbf{S},Y)$.

Since it is still a challenge to optimally choose the bandwidth for finite and especially small sample sizes, if we want to estimate the population mean directly, we rather rely here on the idea of scale spaces, which was originally introduced for regression in Chaudhuri and Marron, (1999) and Chaudhuri and Marron, (2000). The goal is then to show that we can provide SCBs for the population mean simultaneously across different scales. Therefore we define now scale spaces derived from Pristley-Rao smoothing estimators which could, however, be replaced by local polynomial or other appropriate linear smoothers.

Definition 3** (Scale Space Process).**

We define the Scale Space Process with respect to a continuous kernel functions $K:\tilde{\mathcal{S}}\times[h_{0},h_{1}]\rightarrow\mathbb{R}$ with $\tilde{\mathcal{S}}\supset\mathcal{S}$ and $0<h_{0}<h_{1}<\infty$ corresponding to Model (24) as

[TABLE]

with scale mean

[TABLE]

In order to apply the previously presented theory we have to obtain first a functional CLT. The version we present is similar to Theorem 3.2 of Chaudhuri and Marron, (2000) with the difference that we consider the limit with respect to the number of observed curves and include the case of possibly having the (random) measurement points depend on the number of samples, too. The regression version in Chaudhuri and Marron, (2000) only treats the limit of observed measurement points going to infinity.

Theorem 8**.**

Let $Y_{1},...,Y_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}Y$ be a sample from Model (24). Assume further that $Z$ has finite second $\mathcal{C}(\mathcal{S})$-moment, $\mathbf{S}$ is a possibly random $P$ vector, where $P$ is allowed to depend on $N$. Moreover, assume that $\max_{s\in\mathcal{S}}\sigma(s)\leq B<\infty$ and

[TABLE]

exists for all $(s,h),(s^{\prime},h^{\prime})\in\mathcal{S}\times\mathcal{H}$, where the expectations are w.r.t. $\mathbf{S}$. Finally, assume that the smoothing kernel $(s,h)\mapsto K(s,h)$ is $\alpha$-Hölder continuous. Then in $\mathcal{C}(\mathcal{S}\times\mathcal{H})$

[TABLE]

Remark 8**.**

Note that an important example covered by the above Theorem is the case where $\mathbf{S}\subset[0,1]$ is non-random and given by $s_{p}=(p-0.5)/P$ for $p=1,...,P$ for all $n=1,...,N$ and $\varepsilon(s_{1}),...,\varepsilon(s_{P}){\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}\mathcal{N}(0,\eta^{2})$. If $P$ is independent of $N$, the Assumption (25) is trivially satisfied. If $P\rightarrow\infty$ as $N\rightarrow\infty$ it is suffcient that the following integral exists and is finite for all $(s,h),(s^{\prime},h^{\prime})\in\mathcal{S}\times\mathcal{H}$

[TABLE]

in order to have Assumption (25) satisfied, which follows for example if $\mathfrak{c}$ and $\mathfrak{e}$ are continuous.

In order to use Theorem 2 we only have to show that we can estimate the LKCs consistently and the assumptions of the GKF are satisfied. The next Proposition gives sufficient conditions on the kernel $K$ to ensure this.

Proposition 2**.**

Assume the setting of Theorem 8 and additionally assume that the kernel $K\in\mathcal{C}^{3}(\mathcal{S}\times\mathcal{H})$. Define

[TABLE]

and assume that $\mathfrak{r}$ has continuous partial derivatives up to order $3$ and the covariance matrices of

[TABLE]

have rank $2$. Then $\hat{\tilde{\mu}}(s,h)$ and $\hat{\tilde{\sigma}}(s,h)$ satisfy Assumptions (L), (E2) and (G3). Thus, all assumptions of Theorem 2 are satisfied.

Remark 9**.**

Assume the situation of Remark 8. Then the assumption that $\mathfrak{r}$ has continuous partial derivatives up to order $3$ follows directly from the assumption $K\in\mathcal{C}^{3}(\mathcal{S}\times\mathcal{H})$.

7 Simulations

We conduct simulations for evaluation of the performance of different methods for constructing SCBs for the population mean of a signal-plus-noise model. We show that the tGKF is a reliable, fast competitor to bootstrap methods, performing well even for small sample sizes. Additionally, the average width of the tGKF SCBs has less variance in our simulations than bootstrap SCBs. On the other hand, we show that among the bootstrap versions our proposed multiplier-$t$ bootstrap with Rademacher weights performs best for symmetric distributions.

Simulations in this section are always based on $5,000$ Monte Carlo simulations in order to estimate the covering rate of the SCBs. We mostly compare SCBs constructed using the tGKF, the non-parametric bootstrap-$t$ (Boots-t) and the multiplier-t with Gaussian (gMult-t) or Rademacher (rMult-t) multipliers. We always use $5,000$ bootstrap replicates.

7.1 Coverage: Smooth Gaussian Case

This first set of simulations deals with the most favourable case for the tGKF method. We simulate samples from the following smooth signal-plus-noise models (5) (call them Model $A$, $B$, $C$)

[TABLE]

with $\mathbf{a}\sim\mathcal{N}(0,I_{7\times 7})$, $\mathbf{b}\sim\mathcal{N}(0,I_{21\times 21})$ and $\mathbf{c}\sim\mathcal{N}(0,I_{36\times 36})$. Moreover, the vector $\mathbf{K}^{A}(s)$ has entries $K^{A}_{i}(s)=\binom{6}{i}s^{i}(1-s)^{6-i}$ the $(i,6)$-th Bernstein polynomial for $i=0,...,6$, $\mathbf{K}^{B}(s)$ has entries K^{C}_{i}(s)=\exp\big{(}-\frac{\left(s-x_{i}\right)^{2}}{2h_{i}^{2}}\big{)} with $x_{i}=i/21$, $h_{i}=0.04$ for $i<10$, $h_{11}=0.2$ and $h_{i}=0.08$ for $i>10$ and $\mathbf{K}^{C}(s)$ is the vector of all entries from the $6\times 6$-matrix K_{ij}(s)=\exp\big{(}-\frac{\|s-x_{ij}\|^{2}}{2h^{2}}\big{)} with $x_{ij}=(i,j)/6$ with $h=0.06$.

Examples of sample paths of the signal-plus-noise models and the error processes, as well as the simulation results are shown in Figure 2. We simulated samples from Model $A$ and $B$ on an equidistant grid of $200$ points of $[0,1]$. Model $C$ was simulated on an equidistant grid with $50$ points in each dimension.

The simulations suggest that the tGKF method and the rMult-t outperforms all competitors, since even for small sample sizes they achieve the correct covering rates.

7.2 Coverage: Smooth Non-Gaussian Case

The tGKF method is based on a formula valid for Gaussian processes only. However, we have seen in Section 4 that under appropriate conditions the covering is expected to be good asymptotically even if this condition is not fulfilled. For the simulations here we use Model $A$ with the only change that $\mathbf{a}$ has i.i.d. entries of a Student’s $t_{3}/\sqrt{3}$ random variable, which is non-Gaussian, but still symmetric. A second simulation study tackles the question of whether non-Gaussian but close to Gaussian processes still will have reasonable covering. Here we use Model $B$ where $\mathbf{b}$ has i.i.d. entries of $(\chi^{2}_{\nu}-\nu)/\sqrt{2\nu}$ random variables for different parameters of $\nu$. These random variables are non-symmetric, but for $\nu\rightarrow\infty$ they converge to a standard normal.

Figure 4 shows that in the symmetric distribution case of Model $A$ the tGKF has some over coverage for small sample sizes, but eventually it is in the targeted region for $95\%$ coverage. We observed that covering is usually quite close to the targeted covering for symmetric non-Gaussian distributions. The case of non-symmetric distributions, however, produces usually undercovering for the tGKF. However, as predicted by Theorem 3 eventually for large $N$ it gets close to the targeted covering. In the case of non-symmetric distributions the bootstrap-t seems to perform better since it converges faster to the correct covering rate and under the symmetry condition the rMult-$t$ still works exceptionally good. This is assumed to be the case since it preserves all moments up to the fourthed.

7.3 Average Width and Variance of Different SCBs

A second important feature of the performance of confidence bands is their average width and the variance of the width, since it is preferable to have the smallest possible width that still has the correct coverage. Moreover, the width of a SCB should remain stable meaning that its variance should be small. For the SCBs discussed in this article the width can be defined as twice the estimated quantile $\hat{q}_{N,\alpha}$ of the corresponding method. Thus, we provide simulations of quantile estimates for various methods for the Gaussian Model $B$ and the non-Gaussian Model $B$ with $\nu=7$. The optimal quantile $q_{\alpha}$ is simulated using a Monte-Carlo simulation with $50,000$ replications.

We compare Degras’ asymptotic method from the R-package SCBmeanfd, the non-parametric bootstrap-$t$, multiplier-$t$, non parametric bootstrap and a simple multiplier bootstrap. Here the latter two methods use the variance of the original sample instead of its bootstrapped version, which we described in Section 3.2.

Table 1 and 2 show the simulation results. We can draw two main conclusions from these tables. First, the tGKF method has the smallest standard error among the compared methods in the widths while still having good coverage (at least asymptotically in the non-Gaussian case). Second, the main competitor –the bootstrap-t– has a much higher standard error in the width, but its coverage converges faster in the highly non-Gaussian and asymmetric case.

7.4 The Influence of Observation Noise

In order to study the influence of observation noise, we evaluate the dependence of the covering rate of the smoothed mean of a signal-plus-noise model with added i.i.d. Gaussian observation noise on the bandwidth and the standard deviation of the observation noise.

For the simulations we generate samples from the Gaussian Model A and Model B on an equidistant grid of $[0,1]$ with 100 points and add $\mathcal{N}(0,\sigma)$-distributed independent observation noise with $\sigma\in\{0.02,0.1,0.2\}$. Afterwards we smooth the samples with a Gaussian kernel with bandwidths $h\in\{0.02,0.03,0.05,0.1\}$. The smoothed curves are evaluated on a equidistant grid with 400 points. The results of these simulations are shown in Figure 5 and 6.

We also study the covering rate of the population mean of the scale process of Model B. The generation of the samples is the same as for the previous simulations. The only difference is that instead of smoothing with one bandwidth we construct the scale space surface. Here we use a equidistant grid of 100 points of the interval $[0.02,0.1]$. The results can be found in Figure 7.

7.5 SCBs for the Difference of Population Means of Two Independent Samples

Since we previously studied the single sample case in great detail, we will only present the case of one dimensional smooth Gaussian noise processes in the two sample scenario. Moreover, we only report the results for the tGKF approach. The previous observations regarding the other methods, however, carry over to this case.

The simulations are designed as follows. We generate two samples of sizes $N$ and $M$ such that $N/M=c\in\{1,2,4\}$. We are interested in four different scenarios. The first scenario is the most favourable having the same correlation structure and the same variance function. Here we use for both samples the Gaussian Model $A$ from subsection 7.1. In all remaining scenarios one of the samples will always be this error process. In order to check the effect of the two samples having different correlation structures, we use Gaussian Model $B$ as the second sample from subsection 7.1. For dependence on the variance, while the correlation structure is the same, we change the variance function in the Gaussian Model $A$ to $\sigma^{2}(s)=0.04$ for the second sample. As an example process where both the correlation and the variance are different we use Gaussian Model $B$ with the modification that the error process has pointwise variance $\sigma^{2}(s)=0.04$. The results of these simulations are shown in Figure 8

8 Application to DTI fibers

We apply our method to a data set studying the impact of the eating disorder anorexia nervosa on white-matter tissue properties during adolescence in young females. The experimental setup and a different methodology to statistically analyze the data using pointwise testing and permutation tests can be found in the original article Travis et al., (2015). The data set consists of $15$ healthy volunteers – a control group – and $15$ patients. For each volunteer $27$ different neural fibers were extracted.

In order to locate differences in the DTI fibers we use two-sample $95\%$-SCBs for the difference in the population mean between the control group and the patients for each fiber over the domain $\mathcal{S}=[0,100]$. Robustness of detected differences across scales is tested by computing also the SCBs for a difference between the population mean of the scale processes as proposed in Section 6. For the latter we used a Gaussian smoothing kernel and the considered bandwidth range $\mathcal{H}=[1.5,10]$ is sampled at $200$ equidistant bandwidths.

The results for the three fibers, where we detect significant differences, are shown in Figure 9. Our results are mostly consistent with the results from (Travis et al., 2015) in the sense that we also detect significant differences in the right thalamic radiation and the left SLF. We also detect differences in the right cingulum hippocampus, which Travis et al did not detect. However, they also found differences for other fibers. Here we have to note that they used several criteria for claiming a significant difference, which did not take simultaneous testing and the correlation structure of this data into account and therefore might be false positives.

Acknowledgments

The authors are very thankful to Alex Bowring and Thomas Nichols for proposing to try Rademacher multipliers instead of Gaussian multipliers, which stimulated some further exploration on the multipliers and improved the performance of the bootstrap method considerably. We are also grateful for the in detail proof reading of Sam Davenport and especially for detection of an error in Theorem $2$.

Appendix A Proofs

A.1 Proof of Proposition 1

Assume that $Z$ is $(\mathcal{L}^{p},\delta)$-Lipschitz, then using convexity of $|\cdot|^{p}$ we compute

[TABLE]

Hence $Z$ has also a finite $p$th $\mathcal{C}(\mathcal{S})$-moment.

A.2 Proof of Claim in Remark 3

Using the multivariate mean value theorem and the Cauchy-Schwarz inequality yields

[TABLE]

Applying the expectation to both sides and then taking the maximum of the resulting sums we obtain

[TABLE]

The proof follows now from the following two observations. Firstly, by Remark 2 each of the expectations we take the maximum of is finite, since all components of the gradient of $\nabla D^{(d,l)}\mathcal{G}$ are Gaussian processes with almost surely continuous sample paths. Secondly, $|\log\|x\||^{-2}\geq x^{2}$ for all $0<x<1$.

A.3 Proof of Theorem 1

Lemma 1**.**

Let the function $\hat{f}_{N}(u)$ and its first derivative $\hat{f}^{\prime}_{N}(u)$ be uniformly consistent estimators of the function $f(u)$ and its first derivative $f^{\prime}(u)$, respectively, where both are uniformly continuous over $u\in\mathbb{R}$. Assume there exists an open interval $I=(a,b)$ such that $f$ is strictly monotone on $I$ and there exists a unique solution $u_{0}\in I$ to the equation $f(u)=0$. Define $\hat{u}_{N}=\sup\{u\in I:\hat{f}_{N}(u)=0\}$. Then $\hat{u}$ is a consistent estimator of $u_{0}$.

Proof.

Assume w.l.o.g. that $f$ is strictly decreasing on $I$. Thus, for any $\varepsilon>0$ we have $f(u_{0}-\varepsilon)>0>f(u_{0}+\varepsilon)$ by $f(u_{0})=0$. The assumption that $\hat{f}_{N}(u)$ is a consistent estimator of $f(u)$ yields

[TABLE]

which implies that with probability tending to 1, there is a root of $\hat{f}$ in $I_{0,\varepsilon}=(u_{0}-\varepsilon,u_{0}+\varepsilon)$. On the other hand the monotonicity of $f$ guarantees the existence of an $\delta>0$ such that $\inf\{|f(u)|:\leavevmode\nobreak\ u\in I\setminus I_{0,\varepsilon}\}>\delta$. Moreover, by the uniform consistency of $\hat{f}$, we have that

[TABLE]

Therefore, using the inequality

[TABLE]

we can conclude that

[TABLE]

which implies that with probability tending to 1, there is no root of $\hat{f}_{N}$ outside $I_{0,\varepsilon}$. Hence from the definition of $\hat{u}_{N}$, it is clear that $\hat{u}_{N}$ is the only root of $\hat{f}_{N}$ in $I$ with probability tending to 1. As an immediate consequence we have that

[TABLE]

which finishes the proof that $\hat{u}$ is a consistent estimator of $u_{0}$. ∎

Lemma 2**.**

Let $\hat{L}_{d}^{N}$ be a consistent estimator of $L_{d}$ and $\mathrm{EEC}_{\mathcal{G}}(u)$ given in equation (12).

$\left\|\mathrm{EEC}_{\mathcal{G}}(u)-\widehat{\mathrm{EEC}}_{t_{N-1}}(u)\right\|_{\infty}\xrightarrow{N\rightarrow\infty}0$* almost surely.* 2. 2.

$\left\|\mathrm{EEC^{\prime}}_{\mathcal{G}}(u)-\widehat{\mathrm{EEC^{\prime}}}_{t_{N-1}}(u)\right\|_{\infty}\xrightarrow{N\rightarrow\infty}0$* almost surely.*

Proof.

Part 1. is a direct consequence of the consistency of the LKC estimates and the observation that the EC densities $\rho^{t_{\nu}}$ of a ${t_{\nu}}$-process with $\nu=N-1$ degrees of freedom converges uniformly to the EC densities of a Gaussian process $\rho^{\mathcal{G}}$ as $N$ tends to infinity, i.e.

[TABLE]

The latter follows from Worsley, (1994, Theorem 5.4), which implies that the uniform convergence of EC densities is implied by the uniform convergence of

[TABLE]

To see this, note that the distance

[TABLE]

fulfills $\lim_{u\rightarrow\pm\infty}h_{\nu}(u)=0$. Thus, there is a $C_{\nu}=\max_{u\in\mathbb{R}}|h_{\nu}(u)|$ by continuity of $h_{\nu}$. Moreover, note that $h_{\nu}(u)\geq h_{\nu+1}(u)$ for $\nu\geq 1$, all $u\in\mathbb{R}$ and $\lim_{\nu\rightarrow\infty}h_{\nu}(u)=0$. Hence, $C_{\nu}$ converges to zero for $\nu\rightarrow\infty$.

Part 2. uses the same ideas. The only difference is that we have to show that the function

[TABLE]

converges to zero uniformly in $u$. ∎

Proof of Theorem 1.

In order to prove the almost sure convergence $\hat{q}_{\alpha,N}\xrightarrow{N\rightarrow\infty}\tilde{q}_{\alpha}$, we only have to note that for $u$ large enough $u\mapsto\mathrm{EEC}_{\mathcal{G}}(u)$ is strictly monotonically decreasing and therefore the combination of Lemma 1 and 2 yields the claim. ∎

A.4 Proof of Theorem 2

Note that by (Taylor et al., 2005, Theorem 4.3) we have for $\mathcal{Z}$ a mean zero Gaussian process over a parameter set $\mathcal{T}$ that

[TABLE]

where $\sigma_{c}^{2}$ is a variance depending on an associated process to $\mathcal{Z}$. This implies that there is a $\tilde{u}$ such that for all $u\geq\tilde{u}$ we have that

[TABLE]

Equipped with this result using the definition $M=\max_{s\in\mathcal{S}}\mathcal{G}(s)$ and $|M|=\max_{s\in\mathcal{S}}|\mathcal{G}(s)|$ we compute

[TABLE]

Here $I\!I$ converges to zero for $N$ tending to infinity by the fCLT for $\hat{\eta}_{N}$ and the consistent estimation of $\hat{\varsigma}_{N}$ from (E1-2). Therefore it remains to treat $I$.

To deal with this summand, note that

[TABLE]

where the Gaussian random process $\mathcal{Z}$ over $\mathcal{T}=\mathcal{S}\times S^{0}$, where $S^{0}=\{1,-1\}$, is defined by $\mathcal{Z}(s,v)=v\cdot\mathcal{G}(s)$.

Using the above equality, $\widehat{\mathrm{EEC}}_{t_{N-1}}(\hat{q}_{\alpha,N})=\alpha/2$, i.e. the definiton of our estimator $\hat{q}_{\alpha,N}$ from Equation (11) and Lemma 2 we have that

[TABLE]

Thus, using the fact that $L_{d}(\mathcal{T},\mathcal{Z})=L_{0}(S^{0})L_{d}(\mathcal{S},\mathcal{G})=2L_{d}(\mathcal{S},\mathcal{G})$ and (27) and the observation that $\tilde{q}_{\alpha}$ is monotonically increasing in $\alpha$ for $\alpha$ small enough, we can bound $I$ by

[TABLE]

for all $\alpha$ smaller than some $\alpha^{\prime}$, which finishes the proof.

∎

Remark 10**.**

The specific definition of $\sigma_{c}$ associated with the Gaussian process $\mathcal{Z}$ can be found in Taylor et al., (2005).

A.5 Proofs of Theorem 3 and 4

The following Lemma provides almost sure uniform convergence results and will be used often in the following proofs.

Lemma 3**.**

Assume that $X$ and $Y$ are $(\mathcal{L}^{1},\delta)$-Lipschitz processes. Let $X_{1},...,X_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}X$ and $Y_{1},...,Y_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}Y$ be two samples. Then $\overline{\textbf{X}}\xrightarrow{N\rightarrow\infty}\mathbb{E}[X]$ uniformly almost surely. If $X$ and $Y$ are $(\mathcal{L}^{2},\delta)$-Lipschitz processes with finite second $\mathcal{C}(\mathcal{S})$-moments. Then $\widehat{\rm cov}_{N}\!\left[\textbf{X},\textbf{Y}\,\right]\xrightarrow{N\rightarrow\infty}{\rm cov}\!\left[X,Y\,\right]$ uniformly almost surely.

Proof.

First claim: Using the generic uniform convergence result Davidson, (1994, Theorem 21.8), we only need to estalish strong stochastical equicontinuity (SSE) of the random function $\overline{\textbf{X}}-\mathbb{E}[X]$, since pointwise convergence is obvious by the SLLNs. SSE, however, can be easily established using (Davidson, 1994, Theorem 21.10 (ii)), since

[TABLE]

for all $s,s^{\prime}\in\mathcal{S}$. Here $A_{1},...,A_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}A$ denote the random variables from the $(\mathcal{L}^{1},\delta)$-Lipschitz property of the $X_{n}$’s and $X$ and hence the random variable $C_{N}$ converges almost surely to the constant $2\mathbb{E}[A]$ by the SLLNs.

Second claim: Adapting the same strategy as above and assuming w.l.o.g. $\mathbb{E}[X]=\mathbb{E}[Y]=0$, we compute

[TABLE]

where $M(X)=\max_{s\in\mathcal{S}}|X(s)|$ and $B_{1},...,B_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}B$ denote the random variables from the $(\mathcal{L}^{2},\delta)$-Lipschitz property of the $Y_{n}$’s and $Y$. Again by the SLLNs the random Lipschitz constant converges almost surely and is finite, since $X$ and $Y$ have finite second $\mathcal{C}(\mathcal{S})$-moments and are $(\mathcal{L}^{2},\delta)-$Lipschitz. ∎

Lemma 4**.**

Let $\mathfrak{c}$ be a covariance function. Then

1. (i)

If $\mathfrak{c}$ is continuous and has continuous partial derivatives up to order $K$, then the Gaussian process $\mathcal{G}(0,\mathfrak{c})$ has $\mathcal{C}^{K}$-sample paths with almost surely uniform and absolutely convergent expansions

[TABLE]

where $\lambda_{i},\varphi_{i}$ are the eigenvalues and eigenfunctions of the covariance operator of $Z$ and $\{A_{i}\}_{i\in\mathbb{N}}$ are i.i.d. $\mathcal{N}(0,1)$. 2. (ii)

*If $Z$ and all its partial derivatives $D^{I}\!Z$ with multi-indices satisfying $|I|\leq K$, $K\in\mathbb{N}$, are *$(\mathcal{L}^{2},\delta)-$Lipschitz processes with finite $\mathcal{C}(\mathcal{S})$-variances, then $\mathfrak{c}$ is continuous and all partial derivatives $D^{I}_{s}\!D^{I^{\prime}}_{s^{\prime}}\!\mathfrak{c}(s,s^{\prime})$ for $I,I^{\prime}\leq K$ exist and are continuous for all $s,s^{\prime}\in\mathcal{S}$.

Proof.

(i) Since $\mathfrak{c}$ is continuous the process $\mathcal{G}(0,\mathfrak{c})$ is mean-square continuous. Hence there is a Karhunen-Loéve expansion of the form

[TABLE]

with $\lambda_{i},\varphi_{i}$ are the eigenvalues and eigenfunctions of the covariance operator associated with $\mathfrak{c}$ and $\{A_{i}\}_{i\in\mathbb{N}}$ are i.i.d. $\mathcal{N}(0,1)$. From Ferreira and Menegatto, (2012, Theorem 5.1) we have that $D^{I}\varphi\in\mathcal{C}^{K}(\mathcal{S})$. Moreover, it is easy to deduce from their equation (4.3) that

[TABLE]

is almost surely absolutely and uniformly convergent. Note that the assumption that $\mathfrak{c}\in\mathcal{C}^{2s}(\mathcal{S}\times\mathcal{S})$ is too strong in their article. They in fact only require in their proofs that all partial derivatives $D^{I}_{s}D^{I^{\prime}}_{s^{\prime}}\mathfrak{c}(s,s^{\prime})$ for $I,I^{\prime}<s$ exist and are continuous.

(ii): The continuity is a simple consequence of the $(\mathcal{L}^{2},\delta)-$Lipschitz property and the finite $\mathcal{C}(\mathcal{S})$-variances. Let $X$ be a process with these properties, then using the Cauchy-Schwarz inequality

[TABLE]

for some $C<\infty$ and therefore $\mathfrak{c}$ and the covariances of $D^{I}\!Z$ are continuous. We only show that $D^{d}_{s}\mathfrak{c}(s,s^{\prime})$ exists and is continuous. The argument is similar for the higher partial derivatives. From the definition we obtain for all $s,s^{\prime}$

[TABLE]

where $e_{d}$ denotes the $d$-th element of the standard basis of $\mathbb{R}^{D}$. Thus, we only have to prove that we can interchange limits and integration. The latter is an immediate consequence of Lebesgue’s dominated convergence theorem, where we obtain the $\mathcal{L}^{1}$ majorant from the $(\mathcal{L}^{2},\delta)-$Lipschitz property as $AZ(s^{\prime})$, where $A\in\mathcal{L}^{2}$. ∎

Proof of Theorem 3.

(i) Since $Z$ is $(\mathcal{L}^{2},\delta)$-Lipschitz the main result in Jain and Marcus, (1975) immediatly implies (E1) with $\tau_{N}=\sqrt{N}$ and $\mathfrak{r}=\mathfrak{c}$. Condition (E2) is obtained from the second part of Lemma 3, since $\sigma(s)Z(s)$ is $(\mathcal{L}^{2},\delta)$-Lipschitz and has finite second $\mathcal{C}(\mathcal{S})$-moment.

(ii) We only need to show that the Gaussian limit process with the covariance $\mathfrak{c}$ fulfills (G1) and (G3). Note that condition (G3) is a consequence of (G1) and the $\mathcal{C}^{3}$-sample paths by Remark 3. But (G1) is already a consequence of Lemma 4. ∎

Lemma 5**.**

Let $Z$ fulfill the assumptions of Theorem 3 (ii) except for (G2) and for all $d,l\in\{1,...,D\}$ suppose that $\mbox{\rm cov}\!\left[\left(D^{d}\!Z(s),D^{(d,l)}\!Z(s)\right)\right]$ has full rank for all $s$. Then $\mathcal{G}=\mathcal{G}(0,\mathfrak{c})$ fullfills (G2).

Proof.

Using the series expansions from Lemma 4 we have that for multi-indices $I_{1},...,I_{P}$, $K\in\mathbb{N}$ and all $v\in\mathbb{R}^{P}$ it follows that

[TABLE]

is convergent for all $s$ (even uniformly). Note that we used here that the expansions are absolutely convergent such that we can change orders in the infinite sums. Thus, it is easy to deduce that $\left(D^{I_{1}}\mathcal{G},...,D^{I_{P}}\mathcal{G}\right)$ is a Gaussian process.

Therefore \big{(}D^{d}\mathcal{G}(s),D^{(d,l)}\mathcal{G}(s)\big{)} is a multivariate Gaussian random variable for all $s\in\mathcal{S}$, which is nondegenerate if and only if its covariance matrix is non-singular. But this is the case by the assumption, since it is identical to the covariance matrix $\mbox{\rm cov}\!\left[\left(D^{d}\!Z(s),D^{(d,l)}\!Z(s)\right)\right]$. ∎

Proof of Theorem 4.

The proof is almost identical to the proof of Theorem 3 and therefore omitted. ∎

A.6 Proof of Theorem 5

First note that using the definition of $R_{n}$ from Equation 18 we obtain

[TABLE]

Thus, the entries of the sample covariance matrix $\hat{\Lambda}_{N}$ are given by

[TABLE]

Now, the second part of Lemma 3 applied to $Z$ and $D^{l}\!Z$ together with the Assumptions (L) and (E2) imply that the first three summands on the r.h.s. converge to zero almost surely in $\mathcal{C}(\mathcal{S})$. Therefore, applying again Lemma 3 to the remaining summand, we obtain

[TABLE]

uniformly almost surely. Thus, $\hat{\Lambda}_{N}\rightarrow\Lambda$ uniformly almost surely.

Assume that $X,Y$ are $(\mathcal{L}^{4},\delta)$-Lipschitz with finite forth $\mathcal{C}(\mathcal{S})$-moments, then

[TABLE]

with $A,B$ the random variables in the $(\mathcal{L}^{4},\delta)$-Lipschitz property of $Y,X$. Note that

[TABLE]

by $(a+b)^{2}\leq 2a^{2}+2b^{2}$ for all $a,b\in\mathbb{R}$ and the Cauchy-Schwarz inequality. Thus, a sample $X_{1}Y_{1},...,X_{N}Y_{N}{\stackrel{{\scriptstyle i.i.d.}}{{\sim}}}XY$ fulfills the assumptions for the CLT in $\mathcal{C}(\mathcal{S})$ given in Jain and Marcus, (1975). Therefore, the following sums converge to a Gaussian process in $\mathcal{C}(\mathcal{S})$:

[TABLE]

Thus, using the latter together with equation (32) and the Assumptions (L) and (E2) we obtain

[TABLE]

This combined with the standard multivariate CLT yields the claim. ∎

A.7 Proof of Theorem 6 and Corollary 1

Proof of Theorem 6.

The almost sure uniform convergence of $\hat{\Lambda}_{N}$ to $\Lambda(s)$ from Theorem 5 implies that the integrands of equations (19), (20) are almost surely uniform convergent, since $\hat{\Lambda}_{N}$ is composed with a differentiable function. Thus, we can interchange the limit $N\rightarrow\infty$ and the integral everywhere except for a set with measure zero. This gives the claim. ∎

Proof of Corollary 1.

We only have to check (L) and (E2) hold true. Therefore note that

[TABLE]

Thus, both uniform convergence results follow from Lemma 3. ∎

A.8 Proof of Theorem 7

We want to use the functional delta method, e.g., Kosorok, (2008, Theorem 2.8). By Theorem 5 the claim follows, if we prove that the corresponding functions are Hadamard differentiable and can compute this derivative.

Case 1D: We have to prove that the function

[TABLE]

is Hadamard differentiable. Therefore, note that the integral is a bounded linear operator and hence it is Fréchet differentiable with derivative being the integral itself. Moreover, $f\mapsto\sqrt{f}$ is Hadamard differentiable by Kosorok, (2008, Lemma 12.2) with Hadamard derivative $DH_{f}(\alpha)=1/\sqrt{4f}\alpha$ tangential even to the Skorohod space $D(\mathcal{S})$. Combining this, we obtain the limit distribution $\sqrt{N}(\hat{L}_{1}^{N}-L_{1})$ from the fCLT for $\hat{\Lambda}_{N}$ given in Theorem 5 to be distributed as

[TABLE]

where $G(s)$ is the asymptotic Gaussian process given in Theorem 5.

Case 2D: The strategy of the proof is the same as in 1D, i.e. we need to calculate the Hadamard (Fréchet) derivative of

[TABLE]

The arguments are the same as before. Thus, using the chain rule and derivatives of matrices with respect to their components the Hadamard derivative evaluated at the process $\mathcal{G}$ is given by

[TABLE]

A.9 Proof of Corollary 2

Note that it is well-known that the covariance function of the derivative of a differentiable process with covariance function $\mathfrak{c}$ is given by $\dot{\mathfrak{c}}(s,s^{\prime})=D^{1}_{s}D^{1}_{s^{\prime}}\mathfrak{c}(s,s^{\prime})$. Moreover, using the moment formula for multivariate Gaussian processes we have that

[TABLE]

Combining this with the observation that the variance of the zero mean Gaussian random variable $\frac{1}{2}\int_{\mathcal{S}}\frac{G(s)}{\sqrt{\Lambda(s)}}ds$ is given by

[TABLE]

yields the claim.

A.10 Proof of Theorem 8

We want to apply Pollard, (1990, Theorem 10.6). Therefore, except for the indices we adapt the notations of that theorem and define the necessary variables. Recall that $\max_{s\in\mathcal{S}}\sigma(s)\leq B<\infty$. We obtain

[TABLE]

We have to establish the assumptions (i), (iii) and (iv) as (v) is trivially satisfied in our case and (ii) is Assumption (25). As discussed in Degras, (2011, p.1759) the manageability (i) follows from the inequality

[TABLE]

if $\|(s,h)-(s^{\prime},h^{\prime})\|<\epsilon^{1/\alpha}$. Assumption (iii) follows since we can compute

[TABLE]

and (iv) is due to

[TABLE]

for all $\epsilon>0$, which follows from the convergence theorem for integrals with monotonically increasing integrands and the fact that by Markov’s inequality

[TABLE]

for fixed $\epsilon>0$.

The weak convergence to a Gaussian process now follows from Pollard, (1990, Theorem 10.6). ∎

A.11 Proof of Proposition 2

The first step is to establish that for each $N$ the process

[TABLE]

which has $\mathcal{C}^{3}$-sample paths, has finite second $\mathcal{C}(\mathcal{S})$-moment with a constant uniformly bounded over all $N$, and the process itself and its first derivatives are $(\mathcal{L}^{2},\delta)$-Lipschitz again uniformly over all $N$, since then the same arguments as in the proof of Corollary 1 and Lemma 3 will yield (L) and (E2) and hence the consistency of the LKC estimation by Theorem 6. Therefore, note that

[TABLE]

This yields using $(a+b)^{2}\leq 2(a^{2}+b^{2})$ that

[TABLE]

where the bound is independent of $N$. Basically, the same argument yields the $(\mathcal{L}^{2},\|\cdot\|)$-Lipschitz property for $\tilde{Z}(s,h)$ and all of its partial derivatives up to order $3$ with a bounding $\mathcal{L}^{2}$ random variable independent of $N$.

The differentiability of the sample paths of $\mathcal{G}(0,\mathfrak{r})$ follows from Lemma 4(i).