# Double bootstrapping for visualising the distribution of descriptive
statistics of functional data

Han Lin Shang111Corresponding address: Department of Actuarial Studies and
Business Analytics, Level 7, 4 Eastern Road, Macquarie University, Sydney, NSW
2109, Australia; Email<EMAIL_ADDRESS>ORCID:
https://orcid.org/0000-0003-1769-6430
Department of Actuarial Studies and Business Analytics
Macquarie University

###### Abstract

We propose a double bootstrap procedure for reducing coverage error in the
confidence intervals of descriptive statistics for independent and identically
distributed functional data. Through a series of Monte Carlo simulations, we
compare the finite sample performance of single and double bootstrap
procedures for estimating the distribution of descriptive statistics for
independent and identically distributed functional data. At the cost of longer
computational time, the double bootstrap with the same bootstrap method
reduces confidence level error and provides improved coverage accuracy than
the single bootstrap. Illustrated by a Canadian weather station data set, the
double bootstrap procedure presents a tool for visualising the distribution of
the descriptive statistics for the functional data.

Keywords: confidence interval calibration; coverage error; iterated bootstrap.

##  1 Introduction

Recent computer technology in data collection and storage allows statisticians
to analyse functional data. The monographs by Ramsay and Silverman (2002,
2005) present state-of-art parametric techniques, while the book by Ferraty
and Vieu (2006) gives a nonparametric treatment for functional data analysis.
For selective reviews on various aspects of functional data analysis, consult
Cuevas (2014), Shang (2014), Morris (2015), Goia and Vieu (2016), Wang et al.
(2016) and Reiss et al. (2017).

In functional data analysis, a key objective is to draw statistical inference
about the distribution of descriptive statistics from realisations generated
from an unknown data generating process at a population level. Although it is
important to have a consistent estimator for descriptive statistics, it is
equally important to estimate the variability associated with the descriptive
statistics, construct confidence intervals (CIs) and carry out hypothesis
tests. When such a problem arises, bootstrapping turns out to be the only
practical alternative (see, e.g., Cuevas et al., 2006; McMurry and Politis,
2011).

Since functional data do not have a well-defined density, nonparametric
bootstrapping is commonly used for independent and identically distributed
(i.i.d.) functional data. In nonparametric bootstrapping, one resamples
$\bm{\mathcal{X}}^{b}=\\{\mathcal{X}_{1}^{b},\dots,\mathcal{X}_{n}^{b}\\}$
drawn randomly with replacement from original functions
$\bm{\mathcal{X}}=\\{\mathcal{X}_{1},\dots,\mathcal{X}_{n}\\}$, so that each
$\mathcal{X}_{i}^{b}$ has probability $n^{-1}$ of being equal to any given one
of the $\mathcal{X}_{j}$’s,

${\bf
P}\left(\mathcal{X}_{i}^{b}=\mathcal{X}_{j}|\bm{\mathcal{X}}\right)=n^{-1},\qquad
1\leq i,j\leq n.$

With a set of bootstrap samples
$\\{\bm{\mathcal{X}}^{1},\cdots,\bm{\mathcal{X}}^{B}\\}$, we can reduce
coverage level error between the nominal and empirical coverage probabilities
by iterating the bootstrap once more, where $B$ denotes the number of
bootstrap samples. Through a nonparametric bootstrap, one resamples
$\bm{\mathcal{X}}^{b\eta}=\\{\mathcal{X}_{1}^{b\eta},\dots,\mathcal{X}_{n}^{b\eta}\\}$
drawn randomly with replacement from bootstrap functions
$\bm{\mathcal{X}}^{b}=\\{\mathcal{X}_{1}^{b},\dots,\mathcal{X}_{n}^{b}\\}$, so
that each $\mathcal{X}_{i}^{b\eta}$ has probability $n^{-1}$ of being equal to
any given one of the $\mathcal{X}_{j}^{b}$’s,

${\bf
P}\big{(}\mathcal{X}_{i}^{b\eta}=\mathcal{X}_{j}^{b}|\bm{\mathcal{X}}^{b},\bm{\mathcal{X}}\big{)}=n^{-1},\qquad
1\leq i,j\leq n.$

We use the single and double bootstrap procedures to study the distribution of
descriptive statistics of functional data denoted by
$\theta=\psi(\mathcal{F}_{0})$, which is a function of the underlying
distribution function $\mathcal{F}_{0}$. The bootstrap principle is to
estimate a function of the population distribution $\mathcal{F}_{0}$ by using
the same function of the bootstrap distribution function based on a random
sample of size $n$ from $\mathcal{F}_{0}$. Following the early work by Martin
(1990a, b), we consider an example of iterated bootstrapping, namely double
bootstrapping, to construct CIs with more accurate coverage probability than
single bootstrapping for estimating the distribution of descriptive statistics
of functional data (see also Shang, 2015).

There is an extensive literature on double bootstrapping for scalar-valued
data (see, e.g., Chang and Hall, 2015). The first article to mention the
double bootstrap was by Hall (1986), followed quickly by Beran (1987, 1988).
Extensive theoretical discussions are given in Hall (1988) and Hall and Martin
(1988). In the statistics literature, double bootstrapping was used to
construct CIs and estimate the underlying distribution (see, e.g., Booth and
Hall, 1994; Booth and Presnell, 1998; Lee and Young, 1999; Chang and Hall,
2015). The double bootstrapping was also incorporated within the errors-in-
variables framework and applied to a calibration problem (see Pešta, 2013). In
the econometrics literature, Davidson and MacKinnon (2002) uses the double
bootstrapping to improve the reliability of bootstrap tests of non-nested
linear regression models. The contribution of this paper is to extend the use
of double bootstrapping to functional data analysis. Since functional data are
often observed with measurement error, we aim to provide accurate uncertainty
quantification for the descriptive statistics of functional data.

The outline of this article is as follows: Section 2 sets notations,
definitions, and introduces a bootstrap method. Section 3 presents the
descriptive statistics, the $L_{2}$ distance metrics used for evaluation and
comparison, and the finite sample performance based on a series of simulation
studies. In Section 4, we apply the single and double bootstrap procedures to
a Canadian weather station data set. Section 5 concludes, along with some
ideas on how the double bootstrap procedure presented here can be further
extended. Additional simulation results are shown in the appendices.

##  2 Bootstrapping

###  2.1 Notation

Random functions are assumed to be sampled from a second-order stochastic
process $\mathcal{X}$ in $\mathcal{L}^{2}[0,\tau]$, where
$\mathcal{L}^{2}[0,\tau]$ is the Hilbert space of square-integrable functions
on the bounded interval $[0,\tau]$. The stochastic process $\mathcal{X}$
satisfies the finite variance condition
$\int^{\tau}_{0}\mathcal{X}^{2}(t)dt<\infty$, inner product $\langle
f,g\rangle=\int^{\tau}_{0}f(t)g(t)dt$ for any two functions, $f$ and
$g\in\mathcal{L}^{2}[0,\tau]$ and induced squared norm
$\|\cdot\|^{2}=\langle\cdot,\cdot\rangle$. All random functions are defined on
a common probability space $(\Omega,A,P)$. The notation
$\mathcal{X}\in\mathcal{L}_{\mathcal{H}}^{2}(\Omega,A,P)$ is used to indicate
$\text{E}(\|\mathcal{X}\|^{p})<\infty$ for some $p>0$, where $\mathcal{H}$
denotes the Hilbert space. When $p=1$, $\mathcal{X}(t)$ has the mean curve
$\mu(t)=\text{E}[\mathcal{X}(t)]$; when $p=2$, a non-negative definite
covariance function is given by

$c_{\mathcal{X}}(t,s)=\text{Cov}[\mathcal{X}(s),\mathcal{X}(t)]=\text{E}\left\\{[\mathcal{X}(s)-\mu(s)][\mathcal{X}(t)-\mu(t)]\right\\},$
(2.1)

for all $s,t\in\mathcal{I}$.

The covariance function $c_{\mathcal{X}}(t,s)$ in (2.1) allows the covariance
operator of $\mathcal{X}$, denoted by $\mathcal{K}_{\mathcal{X}}$, to be
defined as

$\mathcal{K}_{\mathcal{X}}(\phi)(s)=\int_{\mathcal{I}}c_{\mathcal{X}}(t,s)\phi(t)dt.$

Via Mercer’s lemma, there exists an orthonormal sequence $(\phi_{k})$ of
continuous function in $\mathcal{L}^{2}[0,\tau]$ and a non-increasing sequence
$\lambda_{k}$ of positive number, such that

$c_{\mathcal{X}}(t,s)=\sum^{\infty}_{k=1}\lambda_{k}\phi_{k}(t)\phi_{k}(s),\qquad
s,t\in\mathcal{I},$

where $(\lambda_{1},\lambda_{2},\dots)$ are eigenvalues, where
$(\phi_{1}(t),\phi_{2}(t),\dots)$ are orthogonal eigenfunctions.

###  2.2 Bootstrap functions

Denote $n$ realisations of a stochastic process $\mathcal{X}$ evaluated at
$t\in\mathcal{I}$ as $\\{\mathcal{X}_{1}(t),\dots,\mathcal{X}_{n}(t)\\}$. When
functional data are i.i.d., we can obtain the bootstrap sample
$\\{\mathcal{X}_{1}^{b}(t),\dots,\mathcal{X}_{n}^{b}(t)\\}$ directly, by
randomly sampling with replacement from the original functions (see also
Shang, 2015). In practice, we can only observe and evaluate $\mathcal{X}$ at
discretised data points $0\leq t_{1}<t_{2}\cdots<t_{T}\leq\tau$, thus the
discretized bootstrap samples are obtained as
$\\{\bm{\mathcal{X}}^{b}(t_{j})=[\mathcal{X}^{b}_{1}(t_{j}),\mathcal{X}^{b}_{2}(t_{j}),\dots,\mathcal{X}_{n}^{b}(t_{j})]^{\top}$;
$j=1,2,\dots,T\\}$.

To avoid the possible appearance of repeated curves in the bootstrap samples,
Cuevas et al. (2006) and Shang (2015) replaced the standard i.i.d. bootstrap
samples by so-called smooth bootstrap samples, which are drawn from a smooth
empirical distribution function. This can be achieved by adding white noise to
the bootstrap sample, expressed as

$\bm{\mathcal{X}}^{0}(t_{j})=\bm{\mathcal{X}}^{b}(t_{j})+\bm{z}(t_{j}),\qquad
j=1,2,\dots,T,$ (2.2)

where
$\\{\bm{z}(t_{j})=[z_{1}(t_{j}),z_{2}(t_{j}),\dots,z_{n}(t_{j})]^{\top}\\}$;
and $[\bm{z}(t_{1}),\bm{z}(t_{2}),\dots,\bm{z}(t_{T})]$ is normally
distributed with mean $\bm{0}$ and covariance matrix
$\beta\bm{\Sigma}_{\bm{\mathcal{X}}}$; and $\bm{\Sigma}_{\bm{\mathcal{X}}}$ is
the covariance matrix of
$[\mathcal{X}(t_{1}),\mathcal{X}(t_{2}),\dots,\mathcal{X}(t_{T})]$ (see also
Shang, 2015). Following Cuevas et al. (2006), we set $\beta=0.05$.

###  2.3 Confidence level estimation

Let $\theta:=\theta(t)=\psi(\mathcal{F}_{0})$ be a statistic whose
distribution is unknown, where $\psi$ is the function that defines the
parameter $\theta$. Let $c_{n}(\alpha)$ denote the largest
$(1-\alpha)^{\text{th}}$ quantile of the distribution of $\theta$. Suppose
$\mathcal{F}$ is the set of all possible values of parameter $\theta$. Then

$\\{\vartheta\in\mathcal{F}:\theta_{\vartheta}\leq c_{n}(\alpha)\\}$

is a confidence set of level $1-\alpha$ for $\theta$.

Although the distribution of a statistic is unknown, we can approximate it by
an empirical distribution of $\mathcal{F}_{1}$. Let $\widehat{c}_{n}(\alpha)$
denote the largest $(1-\alpha)^{\text{th}}$ quantile of the empirical
distribution of $\theta$, and the confidence level will typically converge to
$1-\alpha$ as the sample size $n$ increases (Beran, 1987). However, in a
finite sample, there is a difference in confidence level error. We aim to
study this error via single and double bootstrapping.

Let $\bm{\mathcal{X}}^{b}$ be a bootstrap sample of size $n$ drawn from the
empirical distribution $\mathcal{F}_{1}$, where each component of
$\bm{\mathcal{X}}^{b}$ is conditionally independent, given the original sample
$\bm{\mathcal{X}}$. Denote $\mathcal{F}_{1}^{b}$ as the estimate of
$\mathcal{F}_{0}$ re-computed from the bootstrap sample
$\bm{\mathcal{X}}^{b}$, and let $\bm{\mathcal{X}}^{b\eta}$ be a double
bootstrap sample of size $n$ drawn from the estimated distribution
$\mathcal{F}_{1}^{b}$, where each component of $\bm{\mathcal{X}}^{b\eta}$ is
conditionally independent given $\bm{\mathcal{X}}$ and $\bm{\mathcal{X}}^{b}$
(Beran, 1987).

Let $D(\widehat{\theta},\theta)$ be the distance of a statistic between the
population and sample levels. Let $D(\widehat{\theta}^{b},\widehat{\theta})$
be the distance of a statistic between the original sample and single
bootstrap samples. Let $D(\widehat{\theta}^{b\eta},\widehat{\theta}^{b})$ be
the distance of a statistic between the single and double bootstrap samples.
The empirical coverage probability of $\theta$ obtained from the single
bootstrap is given as

$\frac{\\#\left\\{D(\widehat{\theta}^{b},\widehat{\theta})>D(\widehat{\theta},\theta)\right\\}}{B_{1}},\qquad
b=1,\dots,B_{1},$

where $B_{1}$ is the number of samples in the single bootstrap. In contrast,
the empirical coverage probability of $\theta$ obtained from the double
bootstrap is given as

$\frac{\\#\left\\{D(\widehat{\theta}^{b\eta},\widehat{\theta}^{b})>D(\widehat{\theta},\theta)\right\\}}{B_{1}B_{2}},$

where $B_{2}$ is the number of samples in the second layer of the double
bootstrap. To reduce both computational time and confidence level error,
$B_{2}$ can equal to one (see, e.g., Chang and Hall, 2015). The difference
between the single and double bootstrapping stems from the distance measure
$D(\widehat{\theta}^{b\eta},\widehat{\theta})$, which plays the role of error
correction to achieve better calibration.

##  3 Simulation study

The single and double bootstrap procedures are used to estimate the
distribution of descriptive statistics of functional data. In Section 3.1, we
review several descriptive statistics that characterise i.i.d. functional
data. Section 3.2 introduces a simulated Gaussian process, while the
evaluation metrics of estimation accuracy are given in Section 3.3. Section
3.4 displays the simulation results, where the performance of single and
double bootstrapping are evaluated and compared based on their differences
between the empirical and nominal coverage probabilities.

###  3.1 Descriptive statistics of functional data

We seek an estimator of the functional median, which allows us to rank a set
of functions based on their location depth, i.e., the distance from the
functional median (the deepest curve). This leads to the idea of functional
depth, which has received much attention in the functional data literature
(see e.g., Cuevas et al., 2006, 2007; Gervini, 2012; Lopéz-Pintado and Romo,
2009). We consider two functional depth measures, namely Fraiman and Muniz’s
(2001) depth and Gervini’s (2012) depth based on small ball probability.

Fraiman and Muniz’s (2001) depth is the oldest functional depth measure. For
each $t\in[0,\tau]$, let $\mathcal{F}_{1}$ be the empirical sampling
distribution of $\\{\mathcal{X}_{1}(t),\dots,\mathcal{X}_{n}(t)\\}$ and let
$Z_{i}(t)$ be the univariate depth of function $\mathcal{X}_{i}(t)$, given by

$I_{i}=\int^{\tau}_{0}Z_{i}(t)dt=\int^{\tau}_{0}1-\Big{|}\frac{1}{2}-\mathcal{F}_{1}[\mathcal{X}_{i}(t)]\Big{|}dt,$

and the values of $I_{i}$ provide a way of ranking the curves from inward to
outward (Shang, 2015). The functional median is the deepest curve with the
largest value of $I_{i}$.

Gervini (2012) considered the set of distances between any two functions,
denoted by $d(\mathcal{X}_{i},\mathcal{X}_{j})$. An observation
$\mathcal{X}_{i}$ is an outlier if it is far away from most other functions.
Given a probability $\alpha\in[0,1]$, Gervini (2012) defines the
$\alpha$-radius $r_{i}$ as the distance between $\mathcal{X}_{i}(t)$ and its
$\lceil\alpha n\rceil^{\text{th}}$ closest observations, where $\lceil
x\rceil$ denotes the integer closest to $x$ from above. Customarily, we
consider $\alpha=0.5$. The rank of $r_{i}$ provides a measure of outlyingness
of $\mathcal{X}_{i}(t)$; the smaller the $r_{i}$ is, the more dense the
$\mathcal{X}_{i}(t)$ is (Shang, 2015). The functional median is the curve with
the smallest $r_{i}$.

Apart from the functional median, we also consider sample versions of the
functional mean and functional variance, given by

$\displaystyle\overline{\mathcal{X}}(t)$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\mathcal{X}_{i}(t),$
$\displaystyle\widehat{V}(t)$
$\displaystyle=\frac{1}{n-1}\sum^{n}_{i=1}\left[\mathcal{X}_{i}(t)-\overline{\mathcal{X}}(t)\right]^{2}$

and the $\gamma$-trimmed functional mean, which is the mean function of the
$100(1-\gamma)\%$ deepest curves. The functional trimmed mean is expressed as

$\frac{1}{n-\lceil\gamma n\rceil}\sum^{n-\lceil\gamma
n\rceil}_{i=1}\mathcal{X}_{i}(t),$

where
$\left(\mathcal{X}_{(1)},\mathcal{X}_{(2)},\dots,\mathcal{X}_{(n)}\right)$ are
the ordered sample curves based on their increasing location depth, and
$\gamma\in\Big{[}0,\frac{n-1}{n}\Big{]}$ is the trimming parameter.
Customarily, we consider $\gamma=0.05$.

###  3.2 Simulation setup

A series of Monte Carlo simulation studies is implemented to evaluate and
compare the finite sample performance of the single and double bootstrap
procedures for estimating the distribution of descriptive statistics of i.i.d.
functional data. For comparison, we consider the functional model previously
studied in Cuevas et al. (2006) and Shang (2015): a Gaussian process
$\mathcal{X}(t)$ with mean $m(t)=0.95\times 10t(1-t)+0.05\times 30t(1-t)$,
$\text{Cov}[\mathcal{X}(t_{i}),\mathcal{X}(t_{j})]=\exp(-|t_{i}-t_{j}|/0.3)$,
and $\text{Var}[\mathcal{X}(t)]=1$.

Given that the data generating process is a Gaussian process, we can draw
samples from a multivariate normal distribution with mean
$[m(t_{1}),m(t_{2})$,$\dots$, $m(t_{T})]^{\top}$ and covariance matrix
$\text{Cov}[\mathcal{X}(t_{i}),\mathcal{X}(t_{j})]=\min(t_{i},t_{j})=t_{i}$.
In practice, we evaluate and compare our methods on a common set of grid
points. In our simulation studies, we have taken 101 equally-spaced grid
points between 0 and 1, for two sample sizes $n=100$ and 300.

###  3.3 Simulation evaluation

To evaluate the performance of the bootstrap procedures, we calculate the
bootstrap CIs of the descriptive statistics. Given a set of original functions
$\\{\mathcal{X}_{1},\dots,\mathcal{X}_{n}\\}$, we draw $B_{1}=B_{2}=399$
bootstrap samples for each of the $R=200$ replications; and the same pseudo-
random seed was used for all the methods to ensure the same simulation
randomness (see also Shang, 2015). For each replication, the $100(1-\delta)\%$
bootstrap CIs of a descriptive statistic $\theta$, are defined by calculating
the cut-off value $D(\mathcal{X}_{1},\dots,\mathcal{X}_{n})$, such that
$100(1-\delta)\%$ of the bootstrap repetitions
$(\widehat{\theta}^{b},b=1,\dots,B_{1})$ are within a distance smaller than
$D(\mathcal{X}_{1},\dots,\mathcal{X}_{n})$. In the double bootstrap, the
$100(1-\delta)\%$ bootstrap CIs of a descriptive statistic are defined by
calculating the cut-off value
$D(\mathcal{X}_{1}^{b},\dots,\mathcal{X}_{n}^{b})$, such that
$100(1-\delta)\%$ of the bootstrap repetitions
$(\widehat{\theta}^{b\eta},\eta=1,\dots,B_{2})$ are within a distance smaller
than $D(\mathcal{X}_{1}^{b},\dots,\mathcal{X}_{n}^{b})$.

We calculate the empirical coverage probability that the target function at
the population level lies within the CIs. As a performance measure, $L_{2}$
distance metrics are used for constructing CIs. They are defined as

$\displaystyle\left\|\theta(t)-\widehat{\theta}(t)\right\|_{2}$
$\displaystyle=\left\\{\int_{t\in\mathcal{I}}\left[\theta(t)-\widehat{\theta}(t)\right]^{2}dt\right\\}^{\frac{1}{2}},$
$\displaystyle\left\|\widehat{\theta}(t)-\widehat{\theta}^{b}(t)\right\|_{2}$
$\displaystyle=\left\\{\int_{t\in\mathcal{I}}\left[\widehat{\theta}(t)-\widehat{\theta}^{b}(t)\right]^{2}dt\right\\}^{\frac{1}{2}},$
$\displaystyle\left\|\widehat{\theta}^{b}(t)-\widehat{\theta}^{b\eta}(t)\right\|_{2}$
$\displaystyle=\left\\{\int_{t\in\mathcal{I}}\left[\widehat{\theta}^{b}(t)-\widehat{\theta}^{b\eta}(t)\right]^{2}dt\right\\}^{\frac{1}{2}}.$

In Appendix A, we present additional simulation results using the $L_{\infty}$
distance metric.

###  3.4 Simulation results

In Figure 1, we plot the nominal (from 0.5 to 0.95 in a step of 0.05) and
empirical coverage probabilities for estimating the functional mean of i.i.d.
functional data using single and double bootstrap procedures. As the sample
size increases from $n=100$ to 300, the empirical coverage probability
improves for all the bootstrap procedures, including the single bootstrap
procedure of Shang (2015). This result is not surprising, as the validity of
bootstrapping relies on a moderate or large sample size (see also McMurry and
Politis, 2011). Subject to the same pseudo-random seed, the double bootstrap
procedure outperforms the single bootstrap procedure for most if not all
confidence levels. This result is not surprising, as iterating the bootstrap
principle reduces the dependence between the probability distribution of the
resample and the unknown data generating process (as outlined in Shang, 2015).
Between bootstrapping the original functions and smoothed functions in (2.2),
there is an advantage to bootstrap the original function.

Figure 1: Empirical and nominal coverage probabilities for estimating the
functional mean, based on $B_{1}=B_{2}=399$ repetitions and $R=200$
replications. See (2.2) for the smooth bootstrap.

In Figures 2 and 3, we present the nominal and empirical coverage
probabilities for estimating the functional median based on the Fraiman and
Muniz’s (2001) depth and the $\alpha$-radius depth. A striking feature is that
the single bootstrap procedure noticeably over-estimates the nominal coverage
probability, especially when the confidence level is large. In contrast, the
double bootstrap procedure corrects the coverage error, and it even slightly
under-estimates the nominal coverage probability. Compared to the functional
mean, in general, the functional median is harder to estimate correctly. This
is because the quantity being estimated, such as a quantile, is sensitive to
the discreteness of empirical distribution function (see, e.g., Bickel and
Freedman, 1981; Beran and Srivastava, 1985; Hall and Martin, 1989; Eaton and
Tyler, 1991). However, the comparison result further demonstrates the
usefulness of the double bootstrap procedure, which can still correct coverage
error to achieve better calibration, as measured by the $L_{2}$ metrics.
Between bootstrapping the original and smoothed functions, it is advantageous
to bootstrap smoothed function.

Figure 2: Empirical and nominal coverage probabilities for estimating the
functional median based on the Fraiman and Muniz’s (2001) depth, using
$B_{1}=B_{2}=399$ repetitions and $R=200$ replications.

Figure 3: Empirical and nominal coverage probabilities for estimating the
functional median based on the $\alpha$-radius depth, using $B_{1}=B_{2}=399$
repetitions and $R=200$ replications.

In Figures 4 and 5, we show the nominal and empirical coverage probabilities
for estimating the functional trimmed mean based on the Fraiman and Muniz’s
(2001) depth and the $\alpha$-radius depth. Similar to the functional mean,
the empirical coverage probability improves for all the bootstrap procedures
as the sample size increases from $n=100$ to 300. Subject to the same pseudo-
random seed, the single bootstrap procedure is outperformed by the double
bootstrap procedure for most, if not all, confidence levels. This result again
demonstrates the advantage of the double bootstrap procedure, which can
correct coverage error. Between bootstrapping the original and smoothed
functions, there is an advantage to bootstrap the original function.

Figure 4: Empirical and nominal coverage probabilities for estimating the
functional trimmed mean based on the Fraiman and Muniz’s (2001) depth, using
$B_{1}=B_{2}=399$ repetitions and $R=200$ replications.

Figure 5: Empirical and nominal coverage probabilities for estimating the
functional trimmed mean based on the $\alpha$-radius depth, using
$B_{1}=B_{2}=399$ repetitions and $R=200$ replications.

In Figure 6, we plot the nominal and empirical coverage probabilities for
estimating functional variance. Similar to the functional mean and trimmed
functional mean, the empirical coverage probability improves for all the
bootstrap procedures as the sample size increases from $n=100$ to 300. When
$n=100$, the inferior estimation accuracy of the single bootstrap procedure is
magnified, because the variance is a squared distance between
$\mathcal{X}_{i}(t)$ and $\overline{\mathcal{X}}(t)$. When $n=300$, the
estimation accuracy of the single bootstrap procedure improves, but the double
bootstrap procedure still outperforms it. This result again demonstrates the
usefulness of the double bootstrap procedure because of its better
calibration. Between bootstrapping the original and smoothed functions, there
is an advantage to bootstrap original functions.

Figure 6: Empirical and nominal coverage probabilities for estimating the
functional variance, based on $B_{1}=B_{2}=399$ repetitions and $R=200$
replications.

When the descriptive statistic includes the functional mean, it is
advantageous to bootstrap the original functions; when the descriptive
statistic includes the functional median, it is advantageous to bootstrap the
smoothed functions. In Appendix B, we perform a sensitivity analysis of
bootstrap replications and examine whether or not the bootstrap replications
may affect the finite sample estimation of the empirical coverage
probabilities at various confidence levels.

##  4 Application to meteorology data

We consider the classic Canadian weather station data, which is available
publicly at the fda package (Ramsay et al., 2020). This data set has been
widely studied by Ramsay and Silverman (2002, 2005), in the area of
descriptive analysis and regression analysis of functional data.

In Figure 7, we plot the temperature change in degrees Celsius throughout a
year, taken from 35 weather stations across Canada. The 35 weather stations
cover the Atlantic, Pacific, Continental and Arctic climate zones. The
functional curves were interpolated from 365 data points, which measure the
daily mean temperature recorded by a weather station averaged over the period
from 1960 to 1994. Through the rainbow plot of Hyndman and Shang (2010), the
colours correspond to the geographic climates of stations; the red lines show
the weather stations located in the comparably warmer regions, whereas the
purple lines show the weather stations located in the colder regions (Shang,
2015).

Figure 7: Averaged Canadian daily temperatures from 1960 to 1994 observed at
35 weather stations. Each curve shows averaged Canadian daily temperatures at
a weather station, not at a particular year.

Via the single and double bootstrap procedures, we aim to visualise the
distribution of the descriptive statistics of the Canadian weather station
data. As an illustration, we apply the single and double bootstrap procedures
to plot the 95% CIs of the functional mean, functional variance, functional
median, and 5% trimmed functional mean in Figure 8. The sample estimates are
shown in solid black lines. While the 95% empirical CIs are shown in solid red
lines for the single bootstrap, the 95% empirical CIs are shown in dotted blue
lines for the double bootstrap.

Figure 8: Ninety-five percent CIs of the descriptive statistics for the
Canadian weather station data, based on $B_{1}=B_{2}=399$ repetitions. The
sample estimates are shown in solid black lines; their 95% CIs obtained from
the single and double bootstrap procedures are shown in solid red lines and
dotted blue lines, respectively.

##  5 Conclusion

We present a double bootstrap procedure for drawing random samples from a set
of i.i.d. functional data to visualise the distribution of descriptive
statistics. Through a series of Monte Carlo simulations, we show reduced
confidence level error and improved coverage probability in comparison to the
single bootstrap procedure, using the same bootstrap method. This result is
not surprising since as has been pointed out by Chang and Hall (2015) and
Shang (2015), iterating the bootstrap principle reduces the dependence between
the probability distribution of the resample and the unknown data generating
process. As the number of sample curves increases from 100 to 300, the
estimation accuracy improves for both bootstrap procedures, more so for the
double bootstrap. Illustrated by the Canadian weather station data set, the
single and double bootstrap procedures produce similar results for estimating
the distributions of various descriptive statistics, but there are noticeable
differences at the edges of function support for the functional median.

There are a few ways in which the paper can be further extended, and we
briefly outline two. Firstly, bootstrapping functional time series is still in
its infancy, and we intend to extend the double bootstrap procedure to analyse
stationary and weakly dependent functional time series (see, e.g., Zhu and
Politis, 2017; Nyarige, 2016; Shang, 2018; Paparoditis, 2018; Paparoditis and
Shang, 2020, for various single bootstrap procedures). Secondly, we aim to
develop bootstrap procedures that can handle stationary long-memory functional
time series (see, e.g., Li et al., 2020a) and non-stationary long-memory
functional time series (see, e.g., Li et al., 2020b).

## Acknowledgement

The author thanks Dr. Yanrong Yang for many comments and suggestions.

## Appendix A: $L_{\infty}$ distance metrics

Similar to Section 3.3 of the main manuscript, we also consider $L_{\infty}$
distance metrics for constructing CIs. The $L_{\infty}$ distance metrics are
defined as

$\displaystyle\left\|\theta(t)-\widehat{\theta}(t)\right\|_{\infty}$
$\displaystyle=\sup_{t\in\mathcal{I}}\left|\theta(t)-\widehat{\theta}(t)\right|,$
$\displaystyle\left\|\widehat{\theta}(t)-\widehat{\theta}^{b}(t)\right\|_{\infty}$
$\displaystyle=\sup_{t\in\mathcal{I}}\left|\widehat{\theta}(t)-\widehat{\theta}^{b}(t)\right|,$
$\displaystyle\left\|\widehat{\theta}^{b}(t)-\widehat{\theta}^{b\eta}(t)\right\|_{\infty}$
$\displaystyle=\sup_{t\in\mathcal{I}}\left|\widehat{\theta}^{b}(t)-\widehat{\theta}^{b\eta}(t)\right|.$

Using the $L_{\infty}$ distance metrics for constructing CIs, Figures A.1 to
A.6 present the finite sample performance between the two bootstrap
procedures.

Figure A.1: Empirical and nominal coverage probabilities for estimating the
functional mean, based on $B_{1}=B_{2}=399$ repetitions and $R=200$
replications. See (2.2) for the smooth bootstrap.

Figure A.2: Empirical and nominal coverage probabilities for estimating the
functional median based on the Fraiman and Muniz’s (2001) depth, using
$B_{1}=B_{2}=399$ repetitions and $R=200$ replications.

Figure A.3: Empirical and nominal coverage probabilities for estimating the
functional median based on the $\alpha$-radius depth, using $B_{1}=B_{2}=399$
repetitions and $R=200$ replications.

Figure A.4: Empirical and nominal coverage probabilities for estimating the
functional trimmed mean based on the Fraiman and Muniz’s (2001) depth, using
$B_{1}=B_{2}=399$ repetitions and $R=200$ replications.

Figure A.5: Empirical and nominal coverage probabilities for estimating the
functional trimmed mean based on the $\alpha$-radius depth, using
$B_{1}=B_{2}=399$ repetitions and $R=200$ replications.

Figure A.6: Empirical and nominal coverage probabilities for estimating the
functional variance, based on $B_{1}=B_{2}=399$ repetitions and $R=200$
replications.

## Appendix B: Sensitivity analysis of bootstrap replications

In the double bootstrap procedure, the choices of $B_{1}$ and $B_{2}$ may
affect the empirical coverage probability. In Table LABEL:tab:sensitivity, we
study four different choices of $B_{1}$ and $B_{2}$ and find that the
empirical coverage probabilities are fairly similar regardless different
combinations of $B_{1}$ and $B_{2}$.

Table B.1: Empirical and nominal coverage probabilities are shown for estimating various descriptive statistics, based on four combinations of bootstrap replications. $B_{1}$ symbolizes the number of bootstrap samples in the first layer of the double bootstrap, whereas $B_{2}$ symbolizes the number of bootstrap samples in the second layer of the double bootstrap. FM denotes the Fraiman and Muniz’s (2001) depth. |  |  | Nominal coverage probability
---|---|---|---
Statistic |  |  | 0.5 | 0.55 | 0.6 | 0.65 | 0.7 | 0.75 | 0.8 | 0.85 | 0.9 | 0.95
Mean | $B_{1}=399$ |  L2 $(n=100)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.89 | 0.95
| $B_{2}=399$ |  Linf $(n=100)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.90 | 0.95
|  |  L2 $(n=300)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.84 | 0.90 | 0.95
|  |  Linf $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.80 | 0.85 | 0.90 | 0.95
| $B_{1}=399$ |  L2 $(n=100)$ | 0.48 | 0.53 | 0.58 | 0.63 | 0.68 | 0.73 | 0.78 | 0.84 | 0.89 | 0.94
| $B_{2}=99$ |  Linf $(n=100)$ | 0.48 | 0.53 | 0.58 | 0.63 | 0.68 | 0.73 | 0.79 | 0.84 | 0.89 | 0.94
|  |  L2 $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.89 | 0.94
|  |  Linf $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.89 | 0.94
| $B_{1}=99$ |  L2 $(n=100)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.89 | 0.95
| $B_{2}=399$ |  Linf $(n=100)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.89 | 0.95
|  |  L2 $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.89 | 0.95
|  |  Linf $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.90 | 0.95
| $B_{1}=99$ |  L2 $(n=100)$ | 0.49 | 0.53 | 0.58 | 0.63 | 0.68 | 0.73 | 0.78 | 0.83 | 0.88 | 0.94
| $B_{2}=99$ |  Linf $(n=100)$ | 0.48 | 0.53 | 0.58 | 0.63 | 0.68 | 0.73 | 0.78 | 0.84 | 0.89 | 0.94
|  |  L2 $(n=300)$ | 0.48 | 0.53 | 0.58 | 0.64 | 0.69 | 0.74 | 0.79 | 0.83 | 0.89 | 0.94
|  |  Linf $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.63 | 0.68 | 0.73 | 0.78 | 0.84 | 0.89 | 0.94
Median | $B_{1}=399$ |  L2 $(n=100)$ | 0.33 | 0.41 | 0.50 | 0.59 | 0.65 | 0.70 | 0.75 | 0.82 | 0.87 | 0.92
(FM) | $B_{2}=399$ |  Linf $(n=100)$ | 0.36 | 0.43 | 0.51 | 0.59 | 0.65 | 0.70 | 0.75 | 0.81 | 0.86 | 0.91
|  |  L2 $(n=300)$ | 0.40 | 0.45 | 0.52 | 0.61 | 0.66 | 0.70 | 0.75 | 0.80 | 0.86 | 0.91
|  |  Linf $(n=300)$ | 0.40 | 0.45 | 0.51 | 0.59 | 0.64 | 0.68 | 0.72 | 0.77 | 0.83 | 0.88
| $B_{1}=399$ |  L2 $(n=100)$ | 0.34 | 0.41 | 0.49 | 0.57 | 0.64 | 0.70 | 0.75 | 0.81 | 0.87 | 0.92
| $B_{2}=99$ |  Linf $(n=100)$ | 0.36 | 0.43 | 0.50 | 0.58 | 0.64 | 0.69 | 0.74 | 0.80 | 0.85 | 0.90
|  |  L2 $(n=300)$ | 0.40 | 0.45 | 0.52 | 0.60 | 0.66 | 0.70 | 0.75 | 0.80 | 0.86 | 0.91
|  |  Linf $(n=300)$ | 0.40 | 0.45 | 0.51 | 0.58 | 0.64 | 0.68 | 0.72 | 0.77 | 0.83 | 0.88
| $B_{1}=99$ |  L2 $(n=100)$ | 0.36 | 0.44 | 0.53 | 0.62 | 0.68 | 0.72 | 0.77 | 0.83 | 0.89 | 0.93
| $B_{2}=399$ |  Linf $(n=100)$ | 0.38 | 0.45 | 0.53 | 0.61 | 0.66 | 0.71 | 0.76 | 0.81 | 0.86 | 0.91
|  |  L2 $(n=300)$ | 0.43 | 0.48 | 0.55 | 0.64 | 0.69 | 0.73 | 0.76 | 0.82 | 0.87 | 0.92
|  |  Linf $(n=300)$ | 0.42 | 0.47 | 0.53 | 0.61 | 0.66 | 0.70 | 0.74 | 0.79 | 0.85 | 0.90
| $B_{1}=99$ |  L2 $(n=100)$ | 0.33 | 0.40 | 0.49 | 0.57 | 0.64 | 0.69 | 0.75 | 0.81 | 0.87 | 0.92
| $B_{2}=99$ |  Linf $(n=100)$ | 0.36 | 0.42 | 0.50 | 0.57 | 0.64 | 0.69 | 0.74 | 0.80 | 0.85 | 0.90
|  |  L2 $(n=300)$ | 0.39 | 0.45 | 0.51 | 0.58 | 0.65 | 0.69 | 0.74 | 0.79 | 0.85 | 0.91
|  |  Linf $(n=300)$ | 0.40 | 0.45 | 0.51 | 0.57 | 0.63 | 0.67 | 0.71 | 0.76 | 0.82 | 0.87
Median | $B_{1}=399$ |  L2 $(n=100)$ | 0.59 | 0.64 | 0.68 | 0.73 | 0.77 | 0.82 | 0.86 | 0.89 | 0.93 | 0.97
($\alpha$-radius) | $B_{2}=399$ |  Linf $(n=100)$ | 0.72 | 0.77 | 0.80 | 0.84 | 0.87 | 0.90 | 0.93 | 0.95 | 0.97 | 0.99
|  |  L2 $(n=300)$ | 0.56 | 0.61 | 0.66 | 0.70 | 0.75 | 0.79 | 0.84 | 0.88 | 0.92 | 0.96
|  |  Linf $(n=300)$ | 0.68 | 0.73 | 0.77 | 0.81 | 0.85 | 0.88 | 0.91 | 0.94 | 0.96 | 0.99
| $B_{1}=399$ |  L2 $(n=100)$ | 0.58 | 0.63 | 0.68 | 0.72 | 0.76 | 0.81 | 0.85 | 0.89 | 0.92 | 0.96
| $B_{2}=99$ |  Linf $(n=100)$ | 0.71 | 0.76 | 0.80 | 0.83 | 0.87 | 0.90 | 0.92 | 0.95 | 0.97 | 0.99
|  |  L2 $(n=300)$ | 0.56 | 0.60 | 0.65 | 0.70 | 0.74 | 0.79 | 0.83 | 0.87 | 0.91 | 0.95
|  |  Linf $(n=300)$ | 0.67 | 0.72 | 0.76 | 0.80 | 0.84 | 0.87 | 0.91 | 0.93 | 0.96 | 0.98
| $B_{1}=99$ |  L2 $(n=100)$ | 0.59 | 0.64 | 0.68 | 0.73 | 0.77 | 0.81 | 0.86 | 0.89 | 0.93 | 0.96
| $B_{2}=399$ |  Linf $(n=100)$ | 0.72 | 0.76 | 0.80 | 0.84 | 0.87 | 0.90 | 0.93 | 0.95 | 0.97 | 0.99
|  |  L2 $(n=300)$ | 0.56 | 0.61 | 0.66 | 0.70 | 0.75 | 0.79 | 0.84 | 0.88 | 0.92 | 0.96
|  |  Linf $(n=300)$ | 0.68 | 0.72 | 0.76 | 0.81 | 0.84 | 0.88 | 0.91 | 0.94 | 0.97 | 0.99
| $B_{1}=99$ |  L2 $(n=100)$ | 0.58 | 0.63 | 0.68 | 0.72 | 0.77 | 0.81 | 0.85 | 0.89 | 0.92 | 0.96
| $B_{2}=99$ |  Linf $(n=100)$ | 0.71 | 0.76 | 0.80 | 0.83 | 0.87 | 0.90 | 0.93 | 0.95 | 0.97 | 0.99
|  |  L2 $(n=300)$ | 0.55 | 0.60 | 0.65 | 0.69 | 0.74 | 0.78 | 0.83 | 0.87 | 0.91 | 0.95
|  |  Linf $(n=300)$ | 0.67 | 0.71 | 0.76 | 0.80 | 0.84 | 0.87 | 0.90 | 0.93 | 0.96 | 0.98
Trimmed | $B_{1}=399$ |  L2 $(n=100)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.76 | 0.81 | 0.86 | 0.90 | 0.95
mean (FM) | $B_{2}=399$ |  Linf $(n=100)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.71 | 0.76 | 0.81 | 0.86 | 0.91 | 0.96
|  |  L2 $(n=300)$ | 0.51 | 0.56 | 0.61 | 0.66 | 0.71 | 0.76 | 0.81 | 0.86 | 0.90 | 0.95
|  |  Linf $(n=300)$ | 0.51 | 0.56 | 0.61 | 0.66 | 0.71 | 0.76 | 0.81 | 0.85 | 0.90 | 0.95
| $B_{1}=399$ |  L2 $(n=100)$ | 0.50 | 0.55 | 0.59 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95
| $B_{2}=99$ |  Linf $(n=100)$ | 0.49 | 0.55 | 0.59 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95
|  |  L2 $(n=300)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95
|  |  Linf $(n=300)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.94
| $B_{1}=99$ |  L2 $(n=100)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95
| $B_{2}=399$ |  Linf $(n=100)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.81 | 0.86 | 0.91 | 0.95
|  |  L2 $(n=300)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.71 | 0.75 | 0.81 | 0.85 | 0.90 | 0.95
|  |  Linf $(n=300)$ | 0.51 | 0.56 | 0.61 | 0.66 | 0.71 | 0.76 | 0.80 | 0.85 | 0.90 | 0.95
| $B_{1}=99$ |  L2 $(n=100)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.70 | 0.75 | 0.80 | 0.85 | 0.89 | 0.94
| $B_{2}=99$ |  Linf $(n=100)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95
|  |  L2 $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.94
|  |  Linf $(n=300)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.94
Trimmed | $B_{1}=399$ |  L2 $(n=100)$ | 0.51 | 0.56 | 0.61 | 0.66 | 0.71 | 0.76 | 0.81 | 0.86 | 0.91 | 0.95
mean | $B_{2}=399$ |  Linf $(n=100)$ | 0.51 | 0.56 | 0.61 | 0.66 | 0.71 | 0.77 | 0.82 | 0.87 | 0.91 | 0.96
($\alpha$-radius) |  |  L2 $(n=300)$ | 0.51 | 0.56 | 0.61 | 0.66 | 0.71 | 0.76 | 0.81 | 0.86 | 0.90 | 0.95
|  |  Linf $(n=300)$ | 0.51 | 0.56 | 0.61 | 0.66 | 0.71 | 0.76 | 0.81 | 0.86 | 0.91 | 0.95
| $B_{1}=399$ |  L2 $(n=100)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95
| $B_{2}=99$ |  Linf $(n=100)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.71 | 0.76 | 0.81 | 0.86 | 0.90 | 0.95
|  |  L2 $(n=300)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.94
|  |  Linf $(n=300)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95
| $B_{1}=99$ |  L2 $(n=100)$ | 0.50 | 0.55 | 0.61 | 0.66 | 0.71 | 0.76 | 0.81 | 0.86 | 0.90 | 0.95
| $B_{2}=399$ |  Linf $(n=100)$ | 0.50 | 0.56 | 0.61 | 0.66 | 0.71 | 0.76 | 0.82 | 0.86 | 0.91 | 0.96
|  |  L2 $(n=300)$ | 0.50 | 0.55 | 0.61 | 0.66 | 0.71 | 0.76 | 0.80 | 0.85 | 0.90 | 0.95
|  |  Linf $(n=300)$ | 0.51 | 0.56 | 0.61 | 0.66 | 0.71 | 0.76 | 0.81 | 0.86 | 0.91 | 0.95
| $B_{1}=99$ |  L2 $(n=100)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.94
| $B_{2}=99$ |  Linf $(n=100)$ | 0.50 | 0.55 | 0.61 | 0.66 | 0.71 | 0.76 | 0.80 | 0.86 | 0.90 | 0.95
|  |  L2 $(n=300)$ | 0.49 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.84 | 0.89 | 0.94
|  |  Linf $(n=300)$ | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.94
Variance | $B_{1}=399$ |  L2 $(n=100)$ | 0.46 | 0.51 | 0.56 | 0.61 | 0.67 | 0.72 | 0.77 | 0.83 | 0.88 | 0.93
| $B_{2}=399$ |  Linf $(n=100)$ | 0.48 | 0.53 | 0.58 | 0.64 | 0.69 | 0.75 | 0.80 | 0.86 | 0.91 | 0.96
|  |  L2 $(n=300)$ | 0.48 | 0.53 | 0.58 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.89 | 0.94
|  |  Linf $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.70 | 0.75 | 0.80 | 0.86 | 0.91 | 0.96
| $B_{1}=399$ |  L2 $(n=100)$ | 0.45 | 0.50 | 0.55 | 0.61 | 0.66 | 0.71 | 0.76 | 0.82 | 0.87 | 0.92
| $B_{2}=99$ |  Linf $(n=100)$ | 0.47 | 0.52 | 0.58 | 0.63 | 0.68 | 0.74 | 0.79 | 0.85 | 0.90 | 0.95
|  |  L2 $(n=300)$ | 0.48 | 0.53 | 0.58 | 0.63 | 0.68 | 0.73 | 0.78 | 0.83 | 0.88 | 0.94
|  |  Linf $(n=300)$ | 0.48 | 0.53 | 0.58 | 0.64 | 0.69 | 0.74 | 0.79 | 0.85 | 0.90 | 0.95
| $B_{1}=99$ |  L2 $(n=100)$ | 0.46 | 0.51 | 0.56 | 0.61 | 0.67 | 0.72 | 0.77 | 0.83 | 0.88 | 0.93
| $B_{2}=399$ |  Linf $(n=100)$ | 0.48 | 0.53 | 0.58 | 0.64 | 0.69 | 0.75 | 0.80 | 0.86 | 0.91 | 0.96
|  |  L2 $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.64 | 0.69 | 0.74 | 0.79 | 0.84 | 0.89 | 0.95
|  |  Linf $(n=300)$ | 0.49 | 0.54 | 0.59 | 0.65 | 0.70 | 0.75 | 0.80 | 0.86 | 0.91 | 0.96
| $B_{1}=99$ |  L2 $(n=100)$ | 0.45 | 0.50 | 0.55 | 0.60 | 0.66 | 0.71 | 0.76 | 0.82 | 0.87 | 0.93
| $B_{2}=99$ |  Linf $(n=100)$ | 0.47 | 0.52 | 0.58 | 0.63 | 0.69 | 0.74 | 0.79 | 0.85 | 0.90 | 0.95
|  |  L2 $(n=300)$ | 0.48 | 0.53 | 0.58 | 0.63 | 0.68 | 0.73 | 0.78 | 0.83 | 0.89 | 0.94
|  |  Linf $(n=300)$ | 0.48 | 0.53 | 0.59 | 0.64 | 0.69 | 0.74 | 0.80 | 0.85 | 0.90 | 0.95

## References

  * (1)
  * Beran (1987) Beran, R. (1987), ‘Prepivoting to reduce level error in confidence sets’, Biometrika 74(3), 457–468.
  * Beran (1988) Beran, R. (1988), ‘Prepivoting test statistics: A bootstrap view of asymptotic refinements’, Journal of the American Statistical Association: Theory and Methods 83(403), 687–697.
  * Beran and Srivastava (1985) Beran, R. and Srivastava, M. S. (1985), ‘Bootstrap tests and confidence regions for functions of a covariance matrix’, The Annals of Statistics 13(1), 95–115.
  * Bickel and Freedman (1981) Bickel, P. J. and Freedman, D. A. (1981), ‘Some asymptotic theory for the bootstrap’, The Annals of Statistics 9(6), 1196–1217.
  * Booth and Hall (1994) Booth, J. G. and Hall, P. (1994), ‘Monte Carlo approximation and the iterated bootstrap’, Biometrika 81(2), 331–340.
  * Booth and Presnell (1998) Booth, J. G. and Presnell, B. (1998), ‘Allocation of Monte Carlo resources for the iterated bootstrap’, Journal of Computational and Graphical Statistics 7(1), 92–112.
  * Chang and Hall (2015) Chang, J. and Hall, P. (2015), ‘Double-bootstrap methods that use a single double-bootstrap simulation’, Biometrika 102(1), 203–214.
  * Cuevas (2014) Cuevas, A. (2014), ‘A partial overview of the theory of statistics with functional data’, Journal of Statistical Planning and Inference 147, 1–23.
  * Cuevas et al. (2006) Cuevas, A., Febrero, M. and Fraiman, R. (2006), ‘On the use of the bootstrap for estimating functions with functional data’, Computational Statistics & Data Analysis 51(2), 1063–1074.
  * Cuevas et al. (2007) Cuevas, A., Febrero, M. and Fraiman, R. (2007), ‘Robust estimation and classification for functional data via projection-based depth notions’, Computational Statistics 22(3), 481–496.
  * Davidson and MacKinnon (2002) Davidson, R. and MacKinnon, J. G. (2002), ‘Fast double bootstrap tests of nonnested linear regression models’, Econometric Review 21(4), 419–429.
  * Eaton and Tyler (1991) Eaton, M. L. and Tyler, D. E. (1991), ‘On Wielandt’s inequality and its application to the asymptotic distribution of the eigenvalues of a random symmetric matrix’, The Annals of Statistics 19(1), 260–271.
  * Ferraty and Vieu (2006) Ferraty, F. and Vieu, P. (2006), Nonparametric Functional Data Analysis, Springer, New York.
  * Fraiman and Muniz (2001) Fraiman, R. and Muniz, G. (2001), ‘Trimmed mean for functional data’, Test 10(2), 419–440.
  * Gervini (2012) Gervini, D. (2012), ‘Outlier detection and trimmed estimation for general functional data’, Statistica Sinica 22(4), 1639–1660.
  * Goia and Vieu (2016) Goia, A. and Vieu, P. (2016), ‘An introduction to recent advances in high/infinite dimensional statistics’, Journal of Multivariate Analysis 146, 1–6.
  * Hall (1986) Hall, P. (1986), ‘On the bootstrap and confidence intervals’, The Annals of Statistics 14(4), 1431–1452.
  * Hall (1988) Hall, P. (1988), ‘Theoretical comparison of bootstrap confidence intervals’, The Annals of Statistics 16(3), 927–953.
  * Hall and Martin (1988) Hall, P. and Martin, M. A. (1988), ‘On bootstrap resampling and iteration’, Biometrika 75(4), 661–671.
  * Hall and Martin (1989) Hall, P. and Martin, M. A. (1989), ‘A note on the accuracy of bootstrap percentile method confidence intervals for a quantile’, Statistics and Probability Letters 8(3), 197–200.
  * Hyndman and Shang (2010) Hyndman, R. J. and Shang, H. L. (2010), ‘Rainbow plots, bagplots, and boxplots for functional data’, Journal of Computational and Graphical Statistics 19(1), 29–45.
  * Lee and Young (1999) Lee, S. M. S. and Young, G. A. (1999), ‘The effect of Monte Carlo approximation on coverage error of double-bootstrap confidence intervals’, Journal of the Royal Statistical Society: Series B 61(2), 353–366.
  * Li et al. (2020a) Li, D., Robinson, P. M. and Shang, H. L. (2020a), ‘Long-range dependent curve time series’, Journal of the American Statistical Association: Theory and Methods 115(530), 957–971.
  * Li et al. (2020b) Li, D., Robinson, P. M. and Shang, H. L. (2020b), Nonstationary fractionally integrated functional time series, Working paper, University of York.  DOI: 10.13140/RG.2.2.20579.09761.
  * Lopéz-Pintado and Romo (2009) Lopéz-Pintado, S. and Romo, J. (2009), ‘On the concept of depth for functional data’, Journal of the American Statistical Association: Theory and Methods 104(486), 718–734.
  * Martin (1990a) Martin, M. A. (1990a), ‘On bootstrap iteration for coverage correction in confidence intervals’, Journal of the American Statistical Association: Theory and Methods 85(412), 1105–1118.
  * Martin (1990b) Martin, M. A. (1990b), On the double bootstrap, in C. Page and R. LePage, eds, ‘Computing Science and Statistics’, Springer, New York, pp. 73–78.
  * McMurry and Politis (2011) McMurry, T. and Politis, D. N. (2011), Resampling methods for functional data, in F. Ferraty and Y. Romain, eds, ‘The Oxford Handbook of Functional Data Analysis’, Oxford University Press, Oxford, pp. 189–209.
  * Morris (2015) Morris, J. S. (2015), ‘Functional regression’, Annual Review of Statistics and Its Application 2, 321–359.
  * Nyarige (2016) Nyarige, E. G. (2016), The bootstrap for the functional autoregressive model FAR(1), PhD thesis, Technische Universität Kaiserslautern.
https://kluedo.ub.uni-kl.de/frontdoor/index/index/year/2016/docId/4410

  * Paparoditis (2018) Paparoditis, E. (2018), ‘Sieve bootstrap for functional time series’, The Annals of Statistics 46(6B), 3510–3538.
  * Paparoditis and Shang (2020) Paparoditis, E. and Shang, H. L. (2020), Bootstrap prediction bands for functional time series, Technical report, University of Cyprus.
https://arxiv.org/abs/2004.03971

  * Pešta (2013) Pešta, M. (2013), ‘Total least squares and bootstrapping with applications in calibration’, Statistics: A Journal of Theoretical and Applied Statistics 47(5), 966–991.
  * Ramsay et al. (2020) Ramsay, J. O., Wickham, H., Graves, S. and Hooker, G. (2020), fda: Functional Data Analysis.  R package version 5.1.9.
https://CRAN.R-project.org/package=fda

  * Ramsay and Silverman (2002) Ramsay, J. and Silverman, B. (2002), Applied Functional Data Analysis, Springer Series in Statistics, New York.
  * Ramsay and Silverman (2005) Ramsay, J. and Silverman, B. (2005), Functional Data Analysis, 2nd edn, Springer Series in Statistics, New York.
  * Reiss et al. (2017) Reiss, P. T., Goldsmith, J., Shang, H. L. and Ogden, R. T. (2017), ‘Methods for scalar-on-function regression’, International Statistical Review 85(2), 228–249.
  * Shang (2014) Shang, H. L. (2014), ‘A survey of functional principal component analysis’, AStA Advances in Statistical Analysis 98(2), 121–142.
  * Shang (2015) Shang, H. L. (2015), ‘Resampling techniques for estimating the distribution of descriptive statistics of functional data’, Communications in Statistics – Simulation and Computation 44(3), 614–635.
  * Shang (2018) Shang, H. L. (2018), ‘Bootstrap methods for functional time series’, Statistics and Computing 28(1), 1–10.
  * Wang et al. (2016) Wang, J.-L., Chiou, J.-M. and Müller, H.-G. (2016), ‘Functional data analysis’, Annual Review of Statistics and Its Applications 3, 257–295.
  * Zhu and Politis (2017) Zhu, T. and Politis, D. (2017), ‘Kernel estimation of nonparametric functional autoregression models and their bootstrap approximation’, Electronic Journal of Statistics 11(2), 2876–2906.


# Energy levels and extension of the Schrödinger operator.

Y. C. Cantelaube e-mail<EMAIL_ADDRESS>
###### Abstract

Although energy levels are often given by solutions of the radial equation
such that u(0) is non zero, and hence by first-order singular functions which
are not eigenfunctions of H, the latter is always considered as the only
operator that gives energy levels. Vibrational levels of diatomic molecules
are a usual example. We show that the operator which has singular
eigenfunctions, or pseudofunctions, that give energy levels, is the operator
whose action on pseudofunctions amounts to the embedding in the distributions
of R3 of their Hamiltonian in R3/{0}. When its eigenfunctions are regular,
this operator amounts to H. Energy levels, which are given by eigenfunctions
of H when u(0) is zero, are thus given in any case by eigenfunctions of this
operator, which is an extension of the Schrödinger operator, but not of the
Hamiltonian.

U.F.R. de Physique, Université Paris Diderot, Bâtiment Condorcet, 75205 Paris
cedex 13, France

##  1 Introduction. Energy levels and radial equation.

In a central potential $\mathcal{V}$(r) energy levels $\mathcal{E}{}_{n}$ are
eigenvalues of the Hamiltonian that belong to normalized eigenfunctions
$\psi_{n}$ =
[$\mathit{w_{n}}\left(r\right)$/r]$Y_{\ell}^{\mu}\left(\theta,\varphi\right)$
which behave at the origin like rℓ$Y_{\ell}^{\mu}\left(\theta,\varphi\right)$,
and hence eigenvalues of the radial equation that belong to normalized
eigensolutions $\mathit{w_{n}}\left(r\right)$ which behave at the origin like
rℓ+1,

$\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+\frac{\ell(\ell+1)\hbar^{2}}{2mr^{2}}+\mathcal{V}(r)\right]w_{n}\left(r\right)=\mathcal{E}{}_{n}\mathit{w_{n}\left(r\right)}$
(1)

As physical potentials $\mathcal{V}$(r) are not exactly known, one must
substitute theoretical and then approximate potentials V(r) for which the
Hamiltonian has a discrete set of normalizable eigenfunctions $\Psi_{n}$=
[$\mathit{u_{n}}\left(r\right)$/r]$Y_{\ell}^{\mu}\left(\theta,\varphi\right)$,
and hence the radial equation a discrete set of normalizable analytic
eigensolutions $\mathit{u_{n}}\left(r\right)$which behave at the origin like
rℓ+1,

$\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+\frac{\ell(\ell+1)\hbar^{2}}{2mr^{2}}+V\left(r\right)\right]u_{n}\left(r\right)=E_{n}\mathit{u_{n}\left(r\right)}$
(2)

One assumes indeed that if V(r) $\mathcal{\approx V}$(r), $\Psi_{n}$ will be
similar to $\psi_{n}$, $\mathit{u_{n}}\left(r\right)$ to
$\mathit{w_{n}}\left(r\right)$, and $\mathit{E_{n}\approx\mathcal{E}{}_{n}}$.
The solutions of the radial equations, which are written in the form of
series, $\mathit{u_{n}\left(r\right)}$ = $\mathit{r^{\lambda}}$$\sum_{k\geq
0}{\textstyle a_{k}^{n}r^{k}}$, are given in $\textrm{R}^{3}$/{0} by both
roots $\lambda$= $\ell+1$, and $\lambda$= $-\ell$, but in R3 for $\ell$ > 0
only by the root $\lambda$= $\ell+1$ [1]. These solutions are normalizable and
give, when substituted in $\Psi_{n}$, eigenfunctions of H, or solutions of the
Schrödinger equation. For $\ell$ = 0 the radial equations

$\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+\mathcal{V}\left(r\right)\right]w_{n}\left(r\right)=\mathcal{E}{}_{n}\mathit{w_{n}\left(r\right)}$
(3)

$\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+V\left(r\right)\right]u_{n}\left(r\right)=E_{n}\mathit{u_{n}\left(r\right)}$
(4)

have solutions in $\textrm{R}{}^{3}$ given by both roots, that is, $\lambda$ =
1, and $\lambda$= 0. Both roots give normalizable solutions, but only the
solutions given by the root $\lambda$ = 1 give, when substituted in
$\Psi_{n}$, solutions of the Schrödinger equation [1]. Insofar as energy
levels $\mathcal{E}{}_{n}$ are eigenvalues of the Hamiltonian, the solutions
$\mathit{w_{n}}\left(r\right)$ of (3) behave at the origin like r. In order
for $\Psi_{n}$ to be eigenfunctions of H, one must substitute potentials V(r)
for which the following conditions must be satisfied: i) V(r) must be close to
$\mathcal{V}\left(r\right)$, at least in the neighborhood of the minimum; ii)
Eq.(4) must have a discrete set $\mathit{u_{n}}\left(r\right)$ of normalizable
analytic solutions ; iii) these solutions must behave at the origin like r.

Now, explicit analytic solutions of (4) can be derived only for a limited
number of potential energy functions V(r). When we demand that the first two
conditions be satisfied, this number is much more limited. If we moreover
demand that these solutions behave at the origin like r, the problem cannot
necessarily be solved. This is why one substitutes potentials V(r) close to
$\mathcal{V}\left(r\right)$, for which (4) has a discrete set of normalizable
solutions regardless of their value at the origin. One intuitively expects
indeed that if $V\left(r\right)\approx\mathcal{V}\left(r\right)$ and if (4)
has a discrete set of normalizable solutions $\mathit{u_{n}}\left(r\right)$,
moreover demanding that they behave at the origin like r is a supplementary
condition that is not required to obtain eigenvalues $\mathit{E_{n}}$ of (4)
close to the eigenvalues $\mathcal{E}{}_{n}$ of (3). This is shown by the
following example.

##  2 Vibrational levels of diatomic molecules.

The vibrational levels of diatomic molecules are of the form

$\mathcal{E}{}_{n}=-\mathit{V_{m}+(n+\text{\textonehalf)}}\hbar\omega+C_{2}(n+\text{\textonehalf)}^{2}+C_{3}(n+\text{\textonehalf)}^{3}+...\qquad
n=0,1,2...$

In the Born-Oppenheimer approximation the potential energy of interaction
between the two nuclei of a diatomic molecule is a central potential V(r),
attractive at large distances, repulsive at short distances, with a minimum at
r = $\mathit{r_{m}}$. By expanding the potential in powers of r –
$\mathit{r{}_{m}}$ in the neighborhood of the minimum by neglecting the terms
of order $\geq$ 3, it gives

$V\left(r\right)=\text{\textonehalf}m\omega{}^{2}\left(r-r_{m}\right)^{2}-V_{m}\qquad\qquad\omega=[V\,"(r_{m})/m]^{1/2}$
(5)

This parabolic potential is used as a first approximation, which holds in the
neighborhood of the minimum, and then gives the first levels. The radial
equation for $\ell=0$ involving this potential,

$\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+\text{\textonehalf}m\omega{}^{2}\left(r-r_{m}\right)^{2}-V_{m}\right]u_{n}\left(r\right)=E_{n}\mathit{u_{n}\left(r\right)}$
(6)

has the square integrable solutions and the corresponding eigenvalues

$u_{n}\left(r\right)=N_{n}exp\left[-\text{\textonehalf}\beta^{2}(r-r_{m})^{2}\right]H_{n}\left[\beta(r-r_{m})\right]\qquad\beta^{2}=m\omega/\hbar$

$E_{n}=(n+\text{\textonehalf)}\hbar\omega-V_{m}\qquad\qquad\qquad n=0,1,2...$
(7)

where $\mathit{H_{n}}$ are Hermite polynomials of order n and $\mathit{N_{n}}$
normalization constants. Experimental data show that these eigenvalues give
estimates of the first vibrational levels of most molecules. Now, except for
the case where $\beta r_{m}$ coincides with a zero of $\mathit{H_{n}}$,

$u_{n}\left(0\right)=N_{n}exp(-\text{\textonehalf}\beta^{2}r_{m}{}^{2})H_{n}(\beta
r_{m})\neq 0$ (8)

Analytical forms for V(r) closer to the realistic physical potential have been
proposed, the most frequently used is the Morse potential [2]

$V_{M}(r)=V_{m}\\{exp\left[-2a(r-r_{m})\right]-2exp\left[-a(r-r_{m})\right]\\}$
(9)

For $\ell$ = 0 the square-integrable solutions and the eigenvalues of the
radial equation are [2]

$u_{n}(r)=exp\left[-z(r)/2\right]\left[z(r)\right]^{b/2}L_{n+b}^{b}\left[z(r)\right]$

$z=2d\,exp\left[-a(r-r_{m})\right]\qquad
d=(2mV_{m})^{\text{\textonehalf}}(a\hbar)^{-1}\qquad b=2d-1-2n$

$E_{n}=(n+\text{\textonehalf)}\hbar-\left(n+\text{\textonehalf}\right)^{2}(\hbar^{2}\omega^{2}/4V_{m})-V_{m}\qquad\omega=a(2V_{m}/m)^{\text{\textonehalf}}$

where $\mathit{L_{n+b}^{b}}\left(z\right)$ are generalized Laguerre
polynomials. Experimental data show that these eigenvalues give very accurate
values for the vibrational levels of nearly all molecules, but again the
condition $u{}_{n}$(0) = 0 is not satisfied.

If one considers that by substituting the potential (5) “we are left with the
equation of a one-dimensional harmonic oscillator centered at r =
$\mathit{r_{m}}$” [3,$\textrm{A}_{\textrm{V}}$], no boundary condition is
required at the origin, so that in particular perturbation theory for a one-
dimensional harmonic oscillator is applied to the third-order term.

Otherwise, arguments have been proposed to justify that the fact that
$u{}_{n}(0)\neq 0$ is of no significance. According to Pauling et al. the
solutions of the radial equation must be zero at r = 0 and +$\infty$, whereas
the solutions of the harmonic oscillator must be zero at r =
$\text{\textendash}\infty$ and +$\infty$, but because of the rapid decrease in
the harmonic oscillator functions outside the classically permitted region,
“it does not introduce a serious error to consider that the two sets of
boundary conditions are practically equivalent” [4, See also 5,6,7]. As noted
by Schiff, the eigenvalues of the radial equation are those of a linear
harmonic oscillator “if the domain of r is extended to $-\infty$” [8], which
amounts to substituting a one-dimensional problem for a three-dimensional
problem. Or else, one considers that, as (6) is the radial equation for zero
angular momentum states, the exact solutions must be rigorously zero at r = 0,
but the solutions (7) are “practically zero at the origin”
[3,$\textrm{F}{}_{\textrm{VII}}$]. It should be noted that one often considers
indeed that the exact solutions of (6) can be obtained only if the variable is
taken from $-\infty$ to +$\infty$, or if the solutions vanish at the origin.
In fact, (7) are the exact solutions of (6), that is, the solutions in R3
($r\geq$0), as is easily checked.

Similarly, the value of the Morse potential is sometimes considered as so
large at r = 0 that the Morse eigenfunctions are “effectively zero for r < 0”
[9]. In fact they are non zero at r = 0, but according to Morse since “in
every case rΨ will be extremely small … this discrepancy will not affect the
values of the energy levels” [2, See also 6,7].

But if the two boundary conditions, $u_{n}\left(0\right)$ = 0 and
$u{}_{n}\left(0\right)\approx 0$, are practically equivalent, it merely means
that the former condition is not required. First-order singular functions are
not indeed “practically eigenfunctions” of H, they are eigenfunctions of H in
$\textrm{R}^{3}$/{0}, but not in $\textrm{R}^{3}$. In other words, whether or
not u(0) is very small, if it is non zero, Ψn behaves at the origin like 1/r,
and what these examples show is thus that energy levels are not affected
because that they are not given by eigenfunctions of H. By changing the
condition at the origin, we change indeed the operator.

##  3 The energy levels and the operator $\mathit{H_{d}}$.

We must thus determine the operator which has singular eigenfunctions of the
form $\Psi_{n}$ =
[$\mathit{u_{n}}\left(r\right)$/r]$Y_{\ell}^{\mu}\left(\theta,\varphi\right)$
and corresponding eigenvalues $\mathit{E_{n}}$, such that for $\ell$ = 0
$\mathit{u_{n}\left(r\right)}$ and $\mathit{E_{n}}$ are solutions of (4) with
$u_{n}\left(0\right)\neq 0$. Now, when $\Psi_{n}$ and $u_{n}$(r), which are of
$\textrm{C}^{\infty}(\textrm{R}$3/{0}), are singular at the origin, they
define in R3 distributions called pseudofunctions, denoted by
$\textrm{Pf}.\Psi_{n}$ and $\textrm{Pf}.u_{n}$(r) [10], and the latter are not
a solution of (2) [1]. The radial equation which has in R3 singular solutions
$\textrm{Pf}.u_{n}$(r) given by the root $\lambda$ = $-\ell$ is the extension
in R3 of (2), it is the equation [1]

$\textrm{Pf}.\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+\frac{\ell(\ell+1)\hbar^{2}}{2mr^{2}}+V\left(r\right)\right]u_{n}\left(r\right)=E_{n}\textrm{Pf}.\mathit{u_{n}\left(r\right)}$
(10)

Its solutions, given by both roots $\lambda$ = $\ell+1$ and $\lambda$ =
$-\ell$, are either functions $\mathit{u_{n}}\left(r\right)$ which behave at
the origin like $r{}^{\ell+1}$, or pseudofunctions
$\textrm{Pf.}\mathit{u_{n}}\left(r\right)$ which behave at the origin like
$\textrm{Pf.}r^{-\ell}$. (The symbol Pf. is used when we have either the
former or the latter). When the solutions of (10) are given by the roots
$\lambda$ = $\ell+1$ and $\lambda$ = 0 (i.e. $\lambda$ = $-\mathit{\ell}$ for
$\ell$ = 0), the symbol Pf. is useless and (10) is written in the form (2).
Just as $\Psi_{n}$ defines in R3 the pseudofunction
$\textrm{Pf}.\Psi_{n}=\textrm{Pf.}[\mathit{u_{n}}\left(r\right)/r]Y_{\ell}^{\mu}\left(\theta,\varphi\right)$,
its Hamiltonien in R3/{0}, which is the Hamiltonian in the sense of the
functions

$H\Psi_{n}=\frac{1}{r}\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+\frac{\ell(\ell+1)\hbar^{2}}{2mr^{2}}+V(r)\right]\mathit{u_{n}}(r)Y_{\ell}^{\mu}\left(\theta,\varphi\right)\qquad\qquad
r>0$

defines in R3 the distribution

$\textrm{Pf}.H\Psi_{n}=\textrm{Pf}.\left\\{\frac{1}{r}\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+\frac{\ell(\ell+1)\hbar^{2}}{2mr^{2}}+V(r)\right]\mathit{u_{n}}(r)Y_{\ell}^{\mu}\left(\theta,\varphi\right)\right\\}$
(11)

Let $\mathit{H_{d}}$ the operator whose action on a pseudofunction
Pf.$\Psi_{n}$ amounts to the embedding in the distributions of R3 of the
Hamiltonian of $\Psi_{n}$ in R3/{0},

$\mathit{H_{d}}\textrm{Pf}.\Psi_{n}=\textrm{Pf}.H\Psi_{n}$ (12)

If Pf.$\mathit{u_{n}\left(r\right)}$ is a solution of (10), whether it is
given by the root $\lambda$ = $\mathit{\ell}$ \+ 1, or by the root $\lambda$ =
$-\mathit{\ell}$, by substituting (10) in (12), we obtain the eigenvalue
equation of $\mathit{H_{d}}$,

$H_{d}\textrm{Pf.}\Psi_{n}=E_{n}\textrm{Pf.}\Psi_{n}$ (13)

where $\mathit{E_{n}}$ is the eigenvalue of (10) that belongs to
Pf.$\mathit{u_{n}}\left(r\right)$.

The normalizable solutions are given by the roots $\lambda$ =
$\mathit{\ell}$\+ 1 $\forall\ell$, and $\lambda$ = 0, that is, $\lambda\in N$.
As they are less singular at the origin than 1/$r{}^{3}$, the symbol Pf. can
be dropped [10], and (13) is written

$H_{d}\Psi_{n}=E_{n}\Psi_{n}\qquad\qquad\lambda\in N$ (14)

In particular for $\ell$ = 0, $\Psi_{n}$ = $(4\pi){}^{-1/2}[u_{n}(r)/r]$,
whether $u_{n}\left(0\right)=0$ if $\lambda$ = 1, or $u_{n}\left(0\right)\neq
0$ if $\lambda$ = 0, in both cases

$H_{d}\frac{1}{\sqrt{4\pi}}\frac{u_{n}(r)}{r}=\frac{1}{\sqrt{4\pi}}E_{n}\frac{u_{n}(r)}{r}\qquad\qquad\lambda=0$
(15)

where $\mathit{u_{n}\left(r\right)}$ and $\mathit{E_{n}}$ are solutions of the
radial equation (4). $\mathit{H_{d}}$ is thus the required operator.

In order to compare the operators $\mathit{H_{d}}$ and H, let us recall that
the Hamiltonian of
$\textrm{Pf.}\Psi_{n}{\displaystyle{\textstyle=\textrm{Pf}.r^{\lambda-1}\sum_{k\geq
0}a_{k}^{n}r^{k}}\mathit{Y_{\ell}^{\mu}\left(\theta,\varphi\right)}}$ is the
Hamiltonian in the sense of the distributions given by [1]

$H\textrm{Pf}.\Psi_{n}=\textrm{Pf}.\left\\{\frac{1}{r}\left[-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}+\frac{\ell(\ell+1)\hbar^{2}}{2mr^{2}}+V(r)\right]\mathit{u_{n}}(r)Y_{\ell}^{\mu}\left(\theta,\varphi\right)\right\\}-\frac{\hbar^{2}}{2m}Q_{\lambda,\ell}\left(\delta\right)$

$Q_{\lambda,l}\left(\delta\right)=\sum_{k=0}^{-\lambda}a_{k}\chi_{p(k)}B_{\ell,p(k)}C_{p(k)}\mathit{r{}^{\ell}\mathit{\mathit{Y_{\ell}^{\mu}\left(\theta,\varphi\right)}}}\Delta^{p(k)}\delta\qquad
p=-\frac{k+\lambda-\ell}{2}$

$\chi_{p}=\begin{cases}1\;&if\;p\in N\\\ 0&if\;p\notin N\end{cases}\qquad
B_{\ell,p}=\frac{1-2\ell}{4p+1}\qquad
C_{p}=-\frac{(4p+1)\pi^{3/2}}{2^{2p-1}p!\,\Gamma(p+3/2)}$ (16)

where $\Delta^{p}$ is the iterated Laplacian,$\delta$ the Dirac mass, where
the even, resp. odd, terms of the sum are zero , if $\lambda-\ell$ is odd,
resp. even, and where $\mathit{Q_{\lambda,\ell}\left(\delta\right)}$ $\neq$ 0
if and only if there is at least one coefficient $a_{k}\neq 0$ for which k \+
$\lambda-\ell$ is an even negative integer [1]. As

$(H_{d}-H)\textrm{Pf.}\Psi_{n}=\left[\textrm{Pf}.,H\right]\Psi_{n}=\frac{\hbar^{2}}{2m}Q_{\lambda,\ell}\left(\delta\right)$
(17)

the operators $\mathit{H{}_{d}}$ and H differ when they act on pseudofunctions
for which the operators Pf. and H do not commute, or for which
$\mathit{Q_{\lambda,\ell}\left(\delta\right)}$$\neq$ 0\. Considering that
energy levels are given in any case by the operator H amounts to identify the
operators Hd and H, and hence to take no account of the noncommutation of the
operators Pf. and H. As the operators $\mathit{H{}_{d}}$ and H are equivalent
in $\textrm{R}^{3}$/{0}, and in $\textrm{R}^{3}$ when the Laplacian is the
Laplacian in the sense of the functions, it comes from the confusion between
the Laplacians in $\textrm{R}^{3}$/{0} and in $\textrm{R}^{3}$, or/and between
the Laplacians in the sense of the functions and in the sense of the
distributions [11]. By substituting (10) in (16), we obtain

$H\textrm{Pf}.\Psi_{n}=E_{n}\textrm{Pf}.\Psi_{n}-\frac{\hbar^{2}}{2m}Q_{\lambda,\ell}\left(\delta\right)$
(18)

When the solutions of the radial equation are given by the root $\lambda$ =
$\mathit{\ell}$\+ 1, p = – (k \+ 1)/2 $\notin N$ $\forall k$, so that
$Q_{\ell+}{}_{1,\ell}$ = 0, $\Psi_{n}$ behaves at the origin like
$r{}^{\ell}Y_{\ell}^{\mu}\left(\theta,\varphi\right)$, it is an eigenfunction
of H. When the solutions of the radial are given by the root $\lambda$ =
$-\ell$, as $\mathit{a_{o}\neq}$ 0 and p = $\ell\in N$ for k = 0,
$Q{}_{-\ell,\ell}\neq 0$, $\textrm{Pf}.\Psi_{n}$ behaves at the origin like
$\textrm{Pf}.r{}^{-(\ell+1)}Y_{\ell}^{\mu}\left(\theta,\varphi\right)$, it is
not an eigenfunction of H [1]. As $Q{}_{0,0}=-\sqrt{4\pi}u_{n}(0)\delta$, in
the case of normalizable solutions (18) is written

$H\Psi_{n}=E_{n}\Psi_{n}\qquad\qquad\lambda=\ell+1$

$H\,\frac{1}{\sqrt{4\pi}}\frac{u_{n}(r)}{r}=\frac{1}{\sqrt{4\pi}}E_{n}\frac{u_{n}(r)}{r}+\frac{\hbar^{2}\sqrt{\pi}}{m}u_{n}(0)\delta\qquad\qquad\lambda=0$
(19)

Comparison between (14), (15) and (19) shows that energy levels, which are
given by eigenfunctions of H when $u_{n}$(0) = 0, are given in any case by
eigenfunctions of $\mathit{H_{d}}$.

##  4 The Hamiltonian, the Schrödinger operator and the operator
$\mathit{H_{d}}$.

The differential form

$H=-\frac{\hbar^{2}}{2m}\Delta+V$ (20)

is usually called Hamiltonian or Schrödinger operator. Now, a differential
form such as (20) defines an operator only on a given class of functions, so
that the same differential form can define different operators [See e.g. 12].
Any function Ψ, or pseudofunction Pf.Ψ, belongs to the domain $\mathcal{D}$ on
which the Hamiltonian, denoted by H, is defined by (20). On the other hand,
the operator involved in the Schrödinger equation, or Schrödinger operator,
denoted by $\mathit{H_{s}}$, can be defined on the set $\mathcal{D}{}_{s}$ of
the functions which are any superpositions of eigenfunctions of H. As
$\mathcal{D}{}_{s}$ is a subset of $\mathcal{D}$, and as the two operators are
equivalent on $\mathcal{D}{}_{s}$, that is, $H\Psi\equiv H_{s}\Psi$ if
$\Psi\in\mathcal{D}_{s}$, the Schrödinger operator is a restriction of the
Hamiltonian. As long as one confines oneself to solutions of the Schrödinger
equation, these two operators, which have the same eigenfunctions, are
equivalent, so that the distinction is not necessary, and in fact it is not
made. Nevertheless the distinction can be made with respect to the operator
$\mathit{H}_{d}$.

The latter is defined on the same domain $\mathcal{D}$ as H. Moreover (13) and
(18) show that the eigenfunctions of H, or of $\mathit{H{}_{s}}$, form a
subset of the eigenfunctions of $\mathit{H{}_{d}}$, so that if
$\Psi\in\mathcal{D}_{s}$, then $\mathit{H_{d}}\Psi\equiv H_{s}\Psi$. The
operator $\mathit{H_{d}}$ is thus also an extension of the Schrödinger
operator. But the fact that the eigenvalue equation of $\mathit{H{}_{d}}$ is
an extension of the eigenvalue equation of H to the pseudofunctions such that
$\mathit{Q_{\lambda,\ell}\left(\delta\right)}\neq$ 0 does not mean that
$\mathit{H{}_{d}}$ is an extension of H. Eq.(17) shows indeed that they are
different when $\mathit{Q_{\lambda,\ell}\left(\delta\right)}\neq$ 0, or when
the operators Pf. and H do not commute. It follows that the operators
$\mathit{H}$ and $\mathit{H{}_{d}}$ are two different extensions of the
Schrödinger operator.

The condition for the operator $\mathit{H}$ (and then $\mathit{H}_{s}$) to be
self-adjoint is usually written u(0) = 0. However, according to Merzbacher the
condition that the Hamiltonian must be self-adjoint implies that any two
physically admissible eigensolutions of (2) must satisfy the condition [13]

$\textrm{lim }_{r\rightarrow
0}\left(u_{1}^{*}\frac{du_{2}}{dr}-u_{2}\frac{du_{1}^{*}}{dr}\right)=0$ (21)

This condition is satisfied, either for u(0) = 0 and u ’(0)$\neq$ 0, in which
case $\Psi$ is an eigenfunction of $\mathit{H_{d}}$ and of H, or for u ’(0) =
0 and u(0) $\neq$ 0, in which case $\Psi$ is an eigenfunction of
$\mathit{H_{d}}$, but not of H. It follows that (21) is in fact the condition
for $\mathit{H_{d}}$ to be self-adjoint.

The operator $\mathit{H_{d}}$ can also be defined, by using the linearity of
the operator Pf., from the operator $\textrm{Pf}.\Delta$ introduced in [11].
As

$\textrm{Pf}.H=\textrm{Pf}.(-\frac{\hbar^{2}}{2m}\Delta+V)$

by similarly defining the operator $\Delta_{d}$ as

$\Delta_{d}\textrm{Pf}.\Psi=\textrm{Pf}.\Delta\Psi$ (22)

we have

$H_{d}=-\frac{\hbar^{2}}{2m}\Delta_{d}+V$

##  5 A physically not required condition.

The boundary condition $\mathit{u_{n}}$(0) = 0 is definitely stated as a
condition required to obtain energy levels because it is the condition for
$\Psi_{n}$ to be an eigenfunction of the Hamiltonian, but in practice this
condition is not respected. It can be imposed when the solutions are
determined with the help of numerical or approximation methods, but not when
energy levels are determined with the help of analytic solutions.

In the case of diatomic molecules, in the parabolic like in the Morse
potential, the fact that $\mathit{u_{n}}(0)\neq 0$ is regarded as of no
significance owing to the fact that the origin being some distance outside of
the classically allowed region, $\mathit{u_{n}}$(0) is “very small”. As the
parabolic potential (5) is an expansion of the potential in the neighborhood
of the minimum by neglecting the terms of order higher than 2, it remains
close to the physical potential only in the neighborhood of the minimum. The
Morse potential (9) is closer to the physical potential, and in a region more
extended on both sides of the minimum. Comparison between the measured
vibrational levels $\mathcal{E}{}_{n}$ and the eigenvalues $\mathit{E_{n}}$ of
the radial equations involving these potentials shows that the agreement
between these eigenvalues and the energy levels decreases with n, and that the
number of eigenvalues which are close to energy levels is all the greater as
the region in which V(r) $\approx$ $\mathcal{V}$(r) is more extended. It means
that this agreement essentially depends on $\mathit{u_{n}}$(r) in the region
where V(r) is close to $\mathcal{V}$(r), where one therefore expects that
$\mathit{u_{n}}$(r) is close to $\mathit{w_{n}}$(r), the region where nearly
the whole wave function is concentrated. It follows that if energy levels are
not affected by the fact that $\mathit{u_{n}}$(0) $\neq$ 0, it does not come
from the fact that $\mathit{u_{n}}$(0) is “very small”, but from the fact that
the origin is outside of the region where V(r) $\approx\mathcal{V}(r)$, where
$u_{n}\left(r\right)\approx w_{n}\left(r\right)$, and where the wave function
is concentrated, that is, in a region where $\mathit{u_{n}}$(r) has little or
none influence on the fact that $\mathit{E_{n}}$ is close or not to
$\mathcal{E}{}_{n}$.

Besides, in the case of the potential (5) $\mathit{u_{n}}$(r) always vanishes
inside of the classically allowed region. It means in particular that the
condition $\mathit{u_{n}}$(0) = 0 can be satisfied only inside of this region.

Consider a physical potential $\mathcal{V}{}_{o}\left(r\right)$, and a
supplementary potential, or a perturbation, $\mathcal{W}\left(r\right)$, for
which one substitutes the spherical oscillator
$V_{o}\left(r\right)=\text{\textonehalf}m\omega{}^{2}r^{2}$ and the linear
potential W(r) = – Cr, so that the approximation of the physical potential
$\mathcal{V}\left(r\right)$ = $\mathcal{V}{}_{o}\left(r\right)$ \+
$\mathcal{W}\left(r\right)$ is the potential V(r) = $V_{o}\left(r\right)$ \+
W(r) of the form (5), that is,

$V\left(r\right)=\text{\textonehalf}m\omega{}^{2}r^{2}-Cr=\text{\textonehalf}m\omega{}^{2}\left(r-r_{m}\right)^{2}+V_{m}$

$r_{m}=C/m\omega^{2}\qquad\qquad V_{m}=C^{2}/2m\omega^{2}$ (23)

As C and then $\mathit{r{}_{m}}$ have any values, the origin is or not outside
of the classically permitted region. When $\beta\mathit{r_{m}}$ coincides with
a zero of a Hermite polynomial of order N, that is, $\mathit{H_{N}\left(\beta
r_{m}\right)}$ = 0, then $u_{N}\left(0\right)$ = 0, but $u_{n}\left(0\right)$
$\neq$ 0 for n < N. Insofar as the agreement between the eigenvalues of the
radial equation and the energy levels decreases with n, the eigenvalues
$\mathit{E_{n}}$ with n < N that belong to solutions $u_{n}\left(r\right)$
which do not satisfy the condition $u_{n}\left(0\right)$ = 0 must be closer to
an energy level than $\mathit{E_{N}}$ that belongs to a solution
$u_{N}\left(r\right)$ which satisfies this condition. The latter then cannot
be considered as the condition required to obtain accurate, or better
estimates of energy levels.

In the absence of supplementary potential, or of perturbation, the
approximation of the physical potential is the spherical oscillator,
$\mathcal{V}$(r) $\approx$ $V$(r), with

$V\left(r\right)=\text{\textonehalf}m\omega{}^{2}r^{2}$

$u_{n}\left(r\right)=N_{n}exp(-\text{\textonehalf}\beta^{2}r^{2})H_{n}(\beta
r)\qquad E_{n}=(n+\text{\textonehalf)}\hbar\omega$

$u_{n}\left(0\right)=N_{n}H_{n}(0)$

As Hermite polynomials satisfy the condition $H_{2p}\left(0\right)\neq 0$,
$H_{2p+1}\left(0\right)=0$, half the square-integrable solutions of the radial
equation satisfy the condition $\mathit{u_{n}}$(0) = 0, and hence half the
energy levels are given by eigenvalues of H. If the supplementary potential,
or the perturbation, is a 1D potential $\mathcal{W}$(x) for which one
substitutes the linear potential W(x) = – Cx, the potential V(r) =
$V{}_{o}(r)$ + W(x) used as an approximation of the physical potential
$\mathcal{V}$(r) = $\mathcal{V}{}_{o}$(r) + $\mathcal{W}$(x) is still a
spherical oscillator, that is,

$V\left(\rho\right)=\text{\textonehalf}m\omega{}^{2}r^{2}-Cx=\text{\textonehalf}m\omega{}^{2}\rho^{2}-V_{m}$

$\rho=\left[(x-C/m\omega^{2})^{2}+y^{2}+z^{2}\right]\qquad
V_{m}=C^{2}/2m\omega^{2}$

Half the energy levels are still given by eigenvalues of H. But when the
supplementary potential, or the perturbation, is the 3D potential
$\mathcal{W}$(r) for which one substitutes the linear potential W(r) = – Cr,
the potential used as an approximation of the physical potential is the
potential (23), so that the solutions of the radial equation, which are given
by (7), do not satisfy the condition $\mathit{u_{n}}$(0) = 0, and then that
the energy levels are not given by eigenvalues of H, as it is the case for
diatomic molecules. It shows that whether or not energy levels are given by
eigenvalues of H only depends on the potential used as an approximation of the
physical potential – and hence on the latter.

Perturbation methods, in particular variational methods, apply to the
eigenfunctions of H on which the theory is based. By substituting
$\mathit{H_{d}}$ for H these methods are extended to the eigenfunctions of
$\mathit{H_{d}}$ which behave at the origin like 1/r. It is in agreement with
the fact that energy levels are obtained with the help of such functions.
Similarly, situations which are described with the help of functions which
behave at the origin like 1/r, require, to hold in $\textrm{R}{}^{3}$, and not
only in $\textrm{R}{}^{3}$/{0}), that one substitutes $\mathit{H_{d}}$ for H.
This is the case of ingoing, outgoing and standing spherical waves which are
described with the help of Green functions. As the latter behave at the origin
like 1/r, they are eigenfunctions of H in $\textrm{R}{}^{3}$/{0}), but not in
$\textrm{R}{}^{3}$, so that this description does not hold at the origin,
which is physically unsatisfactory. As Green functions are eigenfunctions of
$\mathit{H_{d}}$ in $\textrm{R}{}^{3}$, by substituting this operator, the
description of spherical waves in terms of Green functions holds in the whole
space.

Moreover, wave functions divergent at the origin, but normalizable, are
encountered in the relativistic theory of the hydrogen atom. As the
relativistic generalization of the Schrödinger equation is the Klein-Gordon
equation which governs the wave function in the absence of electromagnetic
field for a relativistic particle with spin zero, if the Laplacian must be
taken in the sense of the distributions in this equation, one will be
similarly led to substitute for $\Delta$ the operator $\Delta_{d}$ defined in
(22).

##  6 Conclusion.

One usually considers, by putting forward different arguments, that the
resolution of the Schrödinger equation must be supplemented with the boundary
condition $\mathit{u_{n}}$(0) = 0, but when this condition is not satisfied,
one admits, by putting forward other arguments, that this is of no
significance. In fact the resolution of the Schrödinger equation need not be
supplemented with any boundary condition [1], but the fact that the condition
$\mathit{u_{n}}$(0) = 0 is not satisfied is significant, since it means,
whether or not $\mathit{u_{n}}$(0) is very small, that energy levels are not
given by solutions of the Schrödinger equation. As energy levels are obtained
by substituting theoretical and then approximate potentials, it should be
taken for granted that in approximate solutions the boundary condition
$\mathit{u_{n}}$(0) = 0, and hence the operator H, or $\mathit{H_{s}}$, are
not required to obtain energy levels, and hence that the operator
$\mathit{H{}_{d}}$, less stringent, is sufficient. The point is that the only
potentials that are exactly known, and then used to obtain energy levels, are
theoretical and then approximate potentials.

## References

  * [1] Y. Cantelaube, Solutions of the Schrödinger equation, boundary condition at the origin, and theory of distributions, arXiv:1203.0551.
  * [2] P. M. Morse, Diatomic molecules according to the wave mechanics. II. Vibrational levels, Phys. Rev. 34, 57 (1929).
  * [3] C. Cohen-Tannoudji, B. Diu and F. Lalöe, Quantum Mechanics (Hermann and John Wiley and sons, Paris-New-York, 1977).
  * [4] L. Pauling and E. Bright Wilson, Introduction to Quantum Mechanics with applications to chemistry (Dover Publications, INC, New-York, 6-415, reprint 1963).
  * [5] H. Eyring, J. Walter and G.E. Kimbali, Quantum Chemistry (John Wiley & Sons, INC., Chapman & Hall, LTD., 1944).
  * [6] I. N. Levine, Molecular spectroscopy, (John Wiley & Sons, 1975).
  * [7] I. N. Levine, Quantum Chemistry, 5th ed. (Prentice Hall, New Jersey, 1991).
  * [8] L. Schiff, Quantum Mechanics, 3rd ed. (McGraw-Hill Book Company, Inc. NY, 1955).
  * [9] K. T. Hecht, Quantum Mechanics (Springer, New York, 2000).
  * [10] L.Schwartz, Théorie des distributions (Paris, Hermann, 1966).
  * [11] Y. Cantelaube and A.L. Khelif, Laplacian in polar coordinates, regular singular function algebra, and theory of distributions, J. Math. Phys. 51, 053518 (2010).
  * [12] E. Prugovecki, Quantum Mechanics in Hilbert Space, 2nd ed. (Academic Press, 1981).
  * [13] E. Merzbacher, Quantum Mechanics, 3rd ed. (John Wiley and sons, INC., 1998).