# Online Estimation for Functional Data

Ying Yang and Fang Yao
Department of Probability and Statistics, School of Mathematical Sciences,
Center for Statistical Science, Peking University, Beijing, China  Fang Yao is
the corresponding author<EMAIL_ADDRESS>This research is supported by
National Natural Science Foundation of China Grants No.11931001 and 11871080,
the LMAM, and the Key Laboratory of Mathematical Economics and Quantitative
Finance (Peking University), Ministry of Education.

###### Abstract

Functional data analysis has attracted considerable interest and is facing new
challenges, one of which is the increasingly available data in a streaming
manner. In this article we develop an online nonparametric method to
dynamically update the estimates of mean and covariance functions for
functional data. The kernel-type estimates can be decomposed into two
sufficient statistics depending on the data-driven bandwidths. We propose to
approximate the future optimal bandwidths by a sequence of dynamically
changing candidates and combine the corresponding statistics across blocks to
form the updated estimation. The proposed online method is easy to compute
based on the stored sufficient statistics and the current data block. We
derive the asymptotic normality and, more importantly, the relative efficiency
lower bounds of the online estimates of mean and covariance functions. This
provides insight into the relationship between estimation accuracy and
computational cost driven by the length of candidate bandwidth sequence.
Simulations and real data examples are provided to support such findings.

Keywords: Dynamic candidate bandwidths; Streaming data; Relative efficiency;
Online learning

##  1 Introduction

Modern technology has promoted the prevalence of functional data in many
fields, such as biomedical studies, engineering, social sciences and so on. A
fundamental problem in functional data analysis is the estimation of mean and
covariance functions, which sets stage for subsequent analyses such as
functional principal component analysis and functional regression.
Technological advancements have improved the speed and volume of data
acquisition and brought new challenges to model such data in a streaming
manner using online approaches, which updates the model estimation with only
the current available data and avoid storing the previous data, as illustrated
in Figure 1. Such online methods not only provide a way to analyze the out-of-
memory data but also allow real-time output of results, and hence are
practically useful. It has been studied extensively in machine learning and
statistics. Langford et al. (2009), Duchi and Singer (2009), Xiao (2010) and
Dekel et al. (2012) discussed stochastic gradient descent problems with
explicit regularizations. Schifano et al. (2016) proposed an online version of
predictive residual test for linear model and improved the aggregated
estimating equation of Lin and Xi (2011) for generalized linear model.
Extensions of linear and quadratic discriminant analyses have also been
studied, see Hiraoka et al. (2000), Kim et al. (2007) and Pang et al. (2005).
Following the custom of these papers, we refer to the classical approaches
using the full data as batch methods.

Figure 1:  Online estimation scheme.

To extend functional data analysis to the online context, the primary interest
is to effectively and efficiently estimate the mean and covariance functions
under various sampling schemes. Suppose that we have $n$ subjects which are
observed at $m_{i}$ time points within the domain $\mathcal{T}$ for the $i$th
individual, denoted by $T_{i}=(T_{i1},\ldots,T_{im_{i}}),i=1,\ldots,n$. In
this work, we focus on the case that each subject is observed at possibly
different times, i.e., $T_{i}$ are randomly distributed, referred to as
independent design. This is distinguished from the common design where all
subjects are observed on a common set of times, as discussed in Cai and Yuan
(2011). From the online perspective, the latter is easier for it may take
sample mean and sample covariance as sufficient statistics and store
historical data in a linear fashion. For independent design, there are
typically two types of estimation methods: pre-smoothing and pooling methods.
The former is to smooth each individual curve and take the cross-sectional
average to form the estimates, which is additive over subjects and can be
extended to online settings straightforwardly. However, pre-smoothing
technique is only applicable to sufficiently dense data and requires $m_{i}\gg
n^{5/4}$, i.e., $m_{i}/n^{5/4}\rightarrow C\in(0,\infty]$, to attain the
root-$n$ convergence rate for the mean and covariance estimation when tuning
parameter (e.g., bandwidth) is chosen optimally for each subject, see Zhang
and Chen (2007) and Kong et al. (2016). By contrast, pooling method can be
adopted to both sparse and dense data. It borrows information from neighboring
data and pools all subjects to estimate the mean and covariance functions,
which is typically achieved by nonparametric methods including splines (Rice
and Silverman, 1991; James et al., 2000; Paul and Peng, 2009) and local
polynomials (Yao et al., 2005a, b; Li and Hsing, 2010), the later only
requires $m_{i}\gg n^{1/4}$ to attain the root-$n$ convergence rate (Zhang and
Wang, 2016).

In this work, we focus on the pooling method using local polynomial smoothing
which is still challenging to extend to online estimation with theoretical
guarantees. Here we take advantage of the extensive theory of local polynomial
regression. Specifically, such kernel-based estimators have the form of
locally weighted least squares and can be decomposed into two sufficient
statistics that are additive on data depending on the bandwidth. The major
difficulty is that the optimal bandwidth changes over data collection and the
sufficient statistics vary accordingly. For instance, Zhou et al. (2003),
Kristan et al. (2010) and Kong and Xia (2019) studied online local polynomial
estimation, but did not update the bandwidths when accessing the stored
statistics from previous data, thus lost the optimality to some extent. To
overcome this, our key proposal is to generate a dynamic sequence of candidate
bandwidths that contain the current optimal bandwidth as well as smaller
values to approximate the future optimal bandwidths. The corresponding
statistics based on the current data block with these dynamically updated
candidate bandwidths are computed, which contain information only from one
data block and are referred to as “sub-statistics”. The stored sufficient
statistics are then updated by aggregating with the current sub-statistics.
Consequently, we are able to update the estimate based on not only the current
data block but also historical data, while only storing one set of dynamically
updated sufficient statistics.

We employ the proposed online mean and covariance estimation to both dense and
sparse functional data, which are also accessible for theoretical analysis. We
derive the asymptotic normality of the online mean and covariance estimates
from sparse to dense designs that is shown to coincide with the phase
transition regime in Zhang and Wang (2016), even when the estimated mean is
used in the covariance estimation. More importantly, when compared to the
classical results based on full data, our online estimates have explicit lower
bounds of the relative efficiency attained by the dynamically updated
bandwidth sequences. These lower bounds indicate that the efficiency is
proportional to the length of candidate bandwidth sequence, which suggests
computing time and memory needed in a linear manner, especially for similar
block sizes. Therefore we can make an informed choice on the trade-off between
estimation accuracy and computational cost according to the real problems at
hand.

The rest of the paper is organized as follows. In Section 2, we briefly
describe the classical estimation of the mean and covariance functions for
functional data analysis. In Section 3, we elucidate the proposed online
method for the mean and covariance estimation using dynamic candidate
bandwidths. In Section 4, we investigate theoretical properties of the
proposed method and establish the lower bounds of the relative efficiency for
the online estimators. Simulations and real data applications are presented in
Section 5 and 6 to validate our theoretical findings. Technical assumptions
and bandwidth selection are delineated in the Appendix, and proofs of
theoretical results are collected in an online Supplementary Material for
space economy.

##  2 Mean and Covariance Estimation

We first review the classical estimation of the mean and covariance functions
for functional data analysis. Let $\\{X(t):t\in\mathcal{T}\\}$ be an $L^{2}$
stochastic process on an interval $\mathcal{T}$ that can be expressed as

$X(t)=\mu(t)+\Phi(t),$

where $\mu(t)={\mathrm{E}}\\{X(t)\\}$ is the mean function, $\Phi(t)$ is the
stochastic part of $X(t)$ with ${\mathrm{E}}\Phi(t)=0$ and $\mathop{\rm
Cov}\\{\Phi(s),\Phi(t)\\}=\gamma(s,t)$ for all $s,t\in\mathcal{T}$, which
implies that $\mathop{\rm Cov}\\{X(s),X(t)\\}=\gamma(s,t)$. In the sequel,
$\mathcal{T}$ is assumed to be $[0,1]$ for simplicity.

Denote the current block as $K$ which increases to $K_{max}$ that tends to
$\infty$ as data accumulate. Suppose that we have observed $n_{k}$ subjects in
the $k$th data block for $k=1,\ldots,K$, among which the $i$th subject has
$m_{ki}$ measurements at possibly irregular time points, $T_{kij}$,
$i=1,\ldots,n_{k}$, $j=1,\ldots,m_{ki}$. In this work, we focus on the setting
where only the number of subjects are increasing as data accumulate. This
setting is different from the online time series analysis, where the
measurements of one or a few sequences are increasing, as discussed in Bacher
et al. (2009), Anava et al. (2013), Yang et al. (2019) and Richard et al.
(2009). They usually considered, e.g., autoregressive models with time-
invariant parameters, hence the accumulation of measurements can improve the
estimation. By contrast, in functional data analysis, the trajectory of each
subject is viewed as a random function, while the mean/covariance functions
are modeled nonparametrically with information accumulated across subjects to
improve estimation. We further remark that, if measurements for each subject
are also collected in a streaming fashion, it raises the challenge due to
dependence between newly arriving measurements and the previously discarded
ones from the same subjects when considering nonparametric modeling of the
mean/covariance function. One possibility is to incorporate some parametric
structures such as those in time series models, coupled with functional data
methods. This leads to a different line of research that can be a topic of
future investigation.

Suppose that observations are contaminated with noises following

$Y_{kij}=X_{ki}(T_{kij})+{\varepsilon}_{kij}=\mu(T_{kij})+\Phi_{ki}(T_{kij})+{\varepsilon}_{kij},$
(1)

where ${\varepsilon}_{kij}$ are independently and identically distributed
(i.i.d) with ${\mathrm{E}}({\varepsilon}_{kij})=0$ and $\mathop{\rm
Var}({\varepsilon}_{kij})=\sigma^{2}$. For notations, the subscript “$k$” or
“$K$” denotes statistics based on only the $k$th or $K$th data block. The
superscript “$(k)$” or “$(K)$” denotes statistics calculated either using the
full data up to $k$ or $K$ by classical batch method (denoted by “
$\widehat{}$ ”), or using only the $k$\- or $K$-th data block and stored
statistics by online method (denoted by “ $\widetilde{}$ ”), both of which
contain information of data blocks up to $k$ or $K$. Recall that Zhang and
Wang (2016) discussed two weighting schemes, equal weight per
observation/subject, which has no material difference when coupled with online
estimation except for their expressions. Thus we focus on the case of equal
weight per observation for conciseness in the sequel.

Denote $W_{h}(\cdot)=W(\cdot/h)/h$, where $W(\cdot)$ is the kernel function
and $h$ is the bandwidth. Classical mean estimate based on local linear
smoother at block $K$ is $\widehat{\mu}^{(K)}(t)=\widehat{\alpha}_{0}$, where

$(\widehat{\alpha}_{0},\widehat{\alpha}_{1})=\mathop{\rm
argmin}\limits_{\alpha_{0},\alpha_{1}}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\sum_{j=1}^{m_{ki}}\\{Y_{kij}-\alpha_{0}-\alpha_{1}(T_{kij}-t)\\}^{2}W_{\widehat{h}_{\mu}^{(K)}}(T_{kij}-t),$

and $\widehat{h}_{\mu}^{(K)}$ is the bandwidth selected by the batch method
using the full data up to block $K$. Let $e_{2}=(1,0)^{\mathrm{T}}$ and
$T_{kij}(t)=(1,T_{kij}-t)^{\mathrm{T}}$, and the solution can be explicitly
written as

$\widehat{\mu}^{(K)}(t)=e_{2}^{\mathrm{T}}\left\\{\sum_{k=1}^{K}P_{k}\big{(}t;\widehat{h}_{\mu}^{(K)}\big{)}\right\\}^{-1}\left\\{\sum_{k=1}^{K}q_{k}\big{(}t;\widehat{h}_{\mu}^{(K)}\big{)}\right\\},$
(2)

where $\\{P_{k},q_{k}\\}$ only depends on the $k$th block given
$\widehat{h}_{\mu}^{(K)}$ by

$\displaystyle
P_{k}(t;\widehat{h}_{\mu}^{(K)})=\sum_{i=1}^{n_{k}}\sum_{j=1}^{m_{ki}}W_{\widehat{h}_{\mu}^{(K)}}(T_{kij}-t)T_{kij}(t)T_{kij}(t)^{\mathrm{T}},$
$\displaystyle
q_{k}(t;\widehat{h}_{\mu}^{(K)})=\sum_{i=1}^{n_{k}}\sum_{j=1}^{m_{ki}}W_{\widehat{h}_{\mu}^{(K)}}(T_{kij}-t)T_{kij}(t)Y_{kij}.$

To estimate $\gamma$, let $m_{ki}\geq 2$ for no contribution to covariance
estimation otherwise. Let $\widehat{C}_{ki,j_{1},j_{2}}^{(K)}$ be the raw
covariance based on $\widehat{\mu}^{(K)}$ in (2), i.e.

$\widehat{C}_{ki,j_{1},j_{2}}^{(K)}=\\{Y_{kij_{1}}-\widehat{\mu}^{(K)}(T_{kij_{1}})\\}\\{Y_{kij_{2}}-\widehat{\mu}^{(K)}(T_{kij_{2}})\\},$
(3)

where $1\leq j_{1}\neq j_{2}\leq m_{ki}$. Let
$\widehat{h}_{\gamma}^{(K)}\in\mathbb{R}$ be the batch bandwidth up to block
$K$ for covariance estimation, and the estimate of $\gamma$ is give by
$\widehat{\gamma}^{(K)}(s,t)=\widehat{\tau}_{0}$, where

$\displaystyle(\widehat{\tau}_{0},\widehat{\tau}_{1},\widehat{\tau}_{2})=\mathop{\rm
argmin}\limits_{\tau_{0},\tau_{1},\tau_{2}}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\qquad\sum_{\mathclap{1\leq
j_{1}\neq j_{2}\leq
m_{ki}}}\big{\\{}\widehat{C}_{ki,j_{1},j_{2}}^{(K)}-\tau_{0}-\tau_{1}(T_{kij_{1}}-s)-\tau_{2}(T_{kij_{1}}-t)\big{\\}}^{2}$
$\displaystyle\times
W_{\widehat{h}_{\gamma}^{(K)}}(T_{kij_{2}}-s)W_{\widehat{h}_{\gamma}^{(K)}}(T_{kij_{2}}-t).$

Let $e_{3}=(1,0,0)^{\mathrm{T}}$ and
$T_{ki,j_{1},j_{2}}(s,t)=(1,T_{kij_{1}}-s,T_{kij_{2}}-t)^{\mathrm{T}}$, and
the explicit solution is

$\widehat{\gamma}^{(K)}(s,t)=e_{3}^{\mathrm{T}}\left\\{\sum_{k=1}^{K}P_{k}\big{(}s,t;\widehat{h}_{\gamma}^{(K)}\big{)}\right\\}^{-1}\left\\{\sum_{k=1}^{K}q_{k}\big{(}s,t;\widehat{h}_{\gamma}^{(K)}\big{)}\right\\},$
(4)

where

$\displaystyle
P_{k}(s,t;\widehat{h}_{\gamma}^{(K)})=\sum_{i=1}^{n_{k}}\qquad\sum_{\mathclap{1\leq
j_{1}\neq j_{2}\leq
m_{ki}}}W_{\widehat{h}_{\gamma}^{(K)}}(T_{kij_{1}}-s)W_{\widehat{h}_{\gamma}^{(K)}}(T_{kij_{2}}-t)T_{ki,j_{1},j_{2}}(s,t)T_{ki,j_{1},j_{2}}(s,t)^{\mathrm{T}},$
$\displaystyle
q_{k}(s,t;\widehat{h}_{\gamma}^{(K)})=\sum_{i=1}^{n_{k}}\qquad\sum_{\mathclap{1\leq
j_{1}\neq j_{2}\leq
m_{ki}}}W_{\widehat{h}_{\gamma}^{(K)}}(T_{kij_{1}}-s)W_{\widehat{h}_{\gamma}^{(K)}}(T_{kij_{2}}-t)T_{ki,j_{1},j_{2}}(s,t)\widehat{C}_{ki,j_{1},j_{2}}.$

##  3 Proposed Online Method

In this section, we give a detailed description of the online estimation for
mean function $\mu(t)$, and the covariance estimation can be derived
similarly. It can be seen from (2) that
$\sum_{k=1}^{K}P_{k}\big{(}t;\widehat{h}_{\mu}^{(K)}\big{)}$ and
$\sum_{k=1}^{K}q_{k}\big{(}t;\widehat{h}_{\mu}^{(K)}\big{)}$ are a pair of
precise sufficient statistics. Given $\widehat{h}_{\mu}^{(K)}$, we only need
to store $\sum_{k=1}^{K}P_{k}\big{(}t;\widehat{h}_{\mu}^{(K)}\big{)}$ and
$\sum_{k=1}^{K}q_{k}\big{(}t;\widehat{h}_{\mu}^{(K)}\big{)}$ in computer
memory instead of the entire $K$ data blocks. However, it is noteworthy that
$\widehat{h}_{\mu}^{(K)}$ takes different values as block $K$ varies, e.g.,
$K\rightarrow\infty$ in streaming problems. This data-driven feature makes the
online local polynomial estimation rather difficult, as storing all
$\big{\\{}\sum_{k=1}^{K}P_{k}\big{(}t;\widehat{h}_{\mu}^{(K)}\big{)},\sum_{k=1}^{K}q_{k}\big{(}t;\widehat{h}_{\mu}^{(K)}\big{)}\big{\\}}_{K=1}^{K_{max}}$
violates the virtue of online method. For example, Kong and Xia (2019) did not
re-calculate the stored statistics with the updated bandwidth. This is a
general issue for online statistical estimation involving data-driven
parameters or structures and partially explains lack of proper online
algorithms for local polynomial estimation that is arguably one of the most
popular nonparametric regression methods.

To overcome this obstacle, we propose to generate a dynamic sequence at block
$k$ ($k\leq K$) consisting of $L$ candidate bandwidths, denoted by
$\\{\eta_{\mu,l}^{(k)}\\}_{l=1}^{L}$. The key idea is to use appropriately
selected $\eta_{\mu,l}^{(k)}$ as a surrogate of future optimal bandwidths
$h_{\mu,*}^{(K)}$ for $k\leq K$. For the flow of presentation, the expressions
of $h_{\mu,*}^{(K)}$ and the corresponding online estimator
$\widetilde{h}_{\mu}^{(k)}$ as well as the optimal form for
$\\{\eta^{(k)}_{\mu,l}\\}_{l=1}^{L}$ are deferred to (17), (18) and (19) of
Section 4. We prove in Theorem 3 that $\widetilde{h}_{\mu}^{(k)}$ converges to
$h_{\mu,*}^{(k)}$ when $k\rightarrow\infty$, and propose to select from
$\\{\eta_{\mu,l}^{(k)}\\}_{l=1}^{L}$ to approximate
$\widetilde{h}_{\mu}^{(K)}$ empirically. Note that $\widetilde{h}_{\mu}^{(k)}$
decreases with respect to $k$, i.e.,
$\widetilde{h}_{\mu}^{(k)}\geq\widetilde{h}_{\mu}^{(K)}$ for $k\leq K$, we set
$\widetilde{h}_{\mu}^{(k)}=\eta_{\mu,1}^{(k)}>\eta_{\mu,2}^{(k)}>\ldots>\eta_{\mu,L}^{(k)}$.
At every block $k$, we compute $L$ pairs of statistics with candidate
bandwidths $\\{\eta_{\mu,l}^{(k)}:l=1,2,\ldots,L\\}$ based on the $k$th data
block. This gives a sequence of statistics
$\\{P_{k}\big{(}t;\eta_{\mu,l}^{(k)}\big{)},q_{k}\big{(}t;\eta_{\mu,l}^{(k)}\big{)}\\}_{l=1}^{L}$
referred to as “sub-statistics” as they contain only information of the $k$th
block (the definition of $P_{k}(t;\cdot),q_{k}(t;\cdot)$ is given in (2)). To
combine information across blocks, we introduce the so-called “pseudo-
sufficient statistics”
$\\{\widetilde{P}_{\mu,l}^{(K)},\widetilde{q}_{\mu,l}^{(K)}\\}_{l=1}^{L}$,
whose update depends on only the sub-statistics of the current data block and
the stored pseudo-sufficient statistics of previous blocks. Specifically,

$\displaystyle\widetilde{P}_{\mu,l}^{(K)}(t)=P_{K}\big{(}t;\eta_{\mu,l}^{(K)}\big{)}+\widetilde{P}_{\mu,{j}_{l}^{(K)}}^{(K-1)}(t),\quad\widetilde{q}_{\mu,l}^{(K)}(t)=q_{K}\big{(}t;\eta_{\mu,l}^{(K)}\big{)}+\widetilde{q}_{\mu,{j}_{l}^{(K)}}^{(K-1)}(t)$
(5)

with initialization $\widetilde{P}_{\mu,l}^{(0)}(t)=0,\
\widetilde{q}_{\mu,l}^{(0)}(t)=0$, where the definition of ${j}_{l}^{(K)}$
would be presented later. The formula (5) indicates that
$\\{\widetilde{P}_{\mu,1}^{(K)},\widetilde{q}_{\mu,1}^{(K)}\\}$ are of form
$\big{\\{}\sum_{k=1}^{K}P_{k}(t;\cdot),\sum_{k=1}^{K}q_{k}(t;\cdot)\big{\\}}$
which are calculated at different candidates across blocks, and hence the
motivation is to select ${j}_{l}^{(K)}$ appropriately to make these candidates
close to $\widetilde{h}_{\mu}^{(K)}$.

Recall that there are $n_{k}$ subjects in the $k$th block and the $i$th one
has $m_{ki}$ measurements. Define

$s_{k,j}=\sum_{i=1}^{n_{k}}m_{ki}(m_{ki}-1)\cdots(m_{ki}-j+1),\
S_{K,j}=\sum_{k=1}^{K}s_{k,j},\ j\geq 1,$ (6)

i.e., $s_{k,1},s_{k,2}$ are the sub-total numbers of observations of the $k$th
block for estimating $\mu$ and $\gamma$, and $S_{K,1},S_{K,2}$ are the total
numbers of observations up to block $K$, respectively. Let
$\omega_{\mu,k}^{(K)}=s_{k,1}/S_{K,1}$ be the weight with respect to numbers
of observations. To update ${j}_{l}^{(K)}$ in an online fashion, we define the
centroids $\\{\phi_{\mu,l}^{(K)}\\}_{l=1}^{L}$ as the weighted average of all
previous candidate bandwidths whose sub-statistics are aggregated, i.e.,

$\phi_{\mu,l}^{(0)}=0,\
\phi_{\mu,l}^{(K)}=(1-\omega_{\mu,K}^{(K)})\phi_{\mu,{j}_{l}^{(K)}}^{(K-1)}+\omega_{\mu,K}^{(K)}\eta_{\mu,l}^{(K)}.$
(7)

We now give the definition of ${j}_{l}^{(K)}$. For $1\leq l\leq L$ and $1\leq
K\leq K_{max}$, we propose to select ${j}_{l}^{(K)}$ to be the index of the
closest one of $\\{\phi_{\mu,l}^{(K-1)}\\}_{l=1}^{L}$ (available at block $K$)
to $\eta_{\mu,l}^{(K)}$, i.e.,

${j}_{l}^{(K)}={\rm
argmin}_{i\in\\{1,2,\ldots,L\\}}|\eta_{\mu,l}^{(K)}-\phi_{\mu,i}^{(K-1)}|.$
(8)

The selection of ${j}_{l}^{(K)}$ in (8) guarantees that
$\phi_{\mu,{j}_{l}^{(K)}}^{(K-1)}$ is close to $\eta_{\mu,l}^{(K)}$, which
makes $\phi_{\mu,l}^{(K)}$ in (7) close to $\eta_{\mu,l}^{(K)}$ as well.
Combined with the definition of $\phi_{\mu,l}^{(K)}$, we conclude that the
combined candidate bandwidths to produce $\phi_{\mu,l}^{(K)}$ are close to
$\eta_{\mu,l}^{(K)}$ on average. Note that
$\eta_{\mu,1}^{(K)}=\widetilde{h}_{\mu}^{(K)}$, the regression estimate at
block $K$ is given by

$\widetilde{\mu}^{(K)}(t)=e_{2}^{\mathrm{T}}\left\\{\widetilde{P}_{\mu,1}^{(K)}(t)\right\\}^{-1}\widetilde{q}_{\mu,1}^{(K)}(t).$
(9)

Since the definitions of ${j}_{l}^{(K)}$ and
$\\{\phi_{\mu,l}^{(K)}\\}_{l=1}^{L}$ are entangled with each other, we take
$K=1,2,8$ and $L=5$ for example in Figure 2 to illustrate the update
procedure. We also display the candidates that produce $\phi_{\mu,i}^{(K)}$
for $i=1,4,7$ with $L=10$ and $K=30,60,100$ in Figure 3. Note that one does
not need to store the previous candidates or the pseudo-bandwidths in
practice.

Figure 2:  A sketch for the update procedure with $L=5$ and $K=1,2,8$. When
$K=1$ (the left panel), $\omega_{\mu,1}^{(1)}=1$ and it can be obtained from
(7) and (5) that $\phi_{\mu,l}^{(1)}=\eta_{\mu,l}^{(1)}$ and
$\widetilde{P}_{\mu,l}^{(1)}(t)=P_{1}\big{(}t;\eta_{\mu,l}^{(1)}\big{)},\widetilde{q}_{\mu,l}^{(1)}(t)=P_{1}\big{(}t;\eta_{\mu,l}^{(1)}\big{)}$
for $1\leq l\leq L$. When $K=2$, $\\{\eta_{\mu,l}^{(2)}\\}_{l=1}^{L}$ are
generated, and for each $l$, ${j}_{l}^{(2)}$ selects the one closest to
$\eta_{\mu,l}^{(2)}$ from $\\{\phi_{\mu,i}^{(1)}\\}_{i=1}^{L}$, then the
corresponding statistics are aggregated as illustrated by the solid lines in
the middle panel. Once ${j}_{l}^{(2)}$ are determined, we update
$\\{\phi_{\mu,l}^{(2)}\\}_{l=1}^{L}$ and
$\big{\\{}\widetilde{P}_{\mu,l}^{(2)},\widetilde{q}_{\mu,l}^{(2)}\big{\\}}_{l=1}^{L}$
according to (7) and (5), and
$\\{{j}_{l}^{(2)},\phi_{\mu,l}^{(1)},\eta_{\mu,l}^{(2)}\\}_{l=1}^{L}$ can be
discarded. For $K=8$, we select ${j}_{l}^{(8)}$ an update
$\\{\phi_{\mu,i}^{(8)}\\}_{i=1}^{L}$ accordingly, and the procedure continues.
By (5),
$\\{\widetilde{P}_{\mu,l}^{(8)}(t),\widetilde{q}_{\mu,l}^{(8)}(t)\\}_{l=1}^{L}$
are comprised of the sub-statistics of blockwise candidates connected by the
dashed lines for the first 8 blocks in the right panel.  Figure 3:  The solid
line by connecting the largest values at each block corresponds to the online
bandwidths $\widetilde{h}_{\mu}^{(k)},k=1,2,\ldots,K$. For illustration, we
set $\widetilde{h}_{\mu}^{(k)}=CS_{k,1}^{-1/5}$ with $C=0.4,S_{k,1}=50k$ for
$k=1,\ldots,100$. The vertical dot series represent the dynamic candidate
sequences $\\{\eta_{\mu,l}^{(k)}\\}_{l=1}^{L}$, $k=1,2,\ldots,K$, which are
calculated by (19) given $\widetilde{h}_{\mu}^{(k)}$. The bold points are the
bandwidths used in
$\big{\\{}\widetilde{P}_{\mu,l}^{(K)},\widetilde{q}_{\mu,l}^{(K)}\big{\\}}$
whose weighted averages are the centroids $\phi_{\mu,l}^{(K)}$ with $l=1,4,7$,
respectively.

To better appreciate (9), let $\widetilde{\eta}_{\mu,k}^{(K)}$ denote the
candidate bandwidth of the $k$th block used in $\widetilde{P}_{\mu,1}^{(K)}$
and $\widetilde{q}_{\mu,1}^{(K)}$, which correspond to the bold dots of $l=1$
in Figure 3 ($K=30,60,100$). We refer to the dynamic sequence
$\\{\widetilde{\eta}_{\mu,k}^{(K)}\\}_{k=1}^{K}$ as the “pseudo-bandwidth”
chain at block $K$ and formulate an equivalent expression of
$\widetilde{\mu}^{(K)}(t)$ in (9):

$\widetilde{\mu}^{(K)}(t)=e_{2}^{\mathrm{T}}\left\\{\sum_{k=1}^{K}P_{k}\big{(}t;\widetilde{\eta}_{\mu,k}^{(K)}\big{)}\right\\}^{-1}\left\\{\sum_{k=1}^{K}q_{k}\big{(}t;\widetilde{\eta}_{\mu,k}^{(K)}\big{)}\right\\},$
(10)

where $P_{k}(t;\cdot),q_{k}(t;\cdot)$ are defined as in (2). Compared (10) to
the batch estimate (2), we are using
$\widetilde{\eta}_{\mu,k}^{(K)},k=1,2,\ldots,K$ to approximate
$\widehat{h}_{\mu}^{(K)}$ which converges to the optimal bandwidth
$h_{\mu,*}^{(K)}$. As mentioned, the definition of centroids and the
combination rule for statistics guarantee that pseudo-bandwidth chain
$\\{\widetilde{\eta}_{\mu,k}^{(K)}\\}_{k=1}^{K}$ is close to
$\eta_{\mu,1}^{(K)}=\widetilde{h}_{\mu}^{(K)}$ on average, as shown in Figure
3. Hence the online estimation using the pseudo-bandwidths shall perform well
as long as $\widetilde{h}_{\mu}^{(K)}$ converges to $h_{\mu,*}^{(K)}$. An
implementation algorithm is given at the end of this section. We mention that
$\\{\eta_{\mu,l}^{(K)}\\}_{l=1}^{L}$ and $\\{{j}_{l}^{(K)}\\}_{l=1}^{L}$ are
recalculated at each $K$, and only the newest
$\\{\phi_{\mu,l}^{(K)}\\}_{l=1}^{L}$ and
$\big{\\{}\widetilde{P}_{\mu,l}^{(K)},\widetilde{q}_{\mu,l}^{(K)}\big{\\}}_{l=1}^{L}$
are stored in memory during the procedure. Hence the algorithm is
computational efficient.

For estimating the covariance function $\gamma(s,t)$, similar to (3), define
the “raw” covariance observations for each block based on the online mean
estimate as

$\widetilde{C}_{ki,j_{1},j_{2}}^{(k)}=\\{Y_{kij_{1}}-\widetilde{\mu}^{(k)}(T_{kij_{1}})\\}\\{Y_{kij_{2}}-\widetilde{\mu}^{(k)}(T_{kij_{2}})\\},$
(11)

where $1\leq j_{1}\neq j_{2}\leq m_{ki},\ i=1,\ldots,n_{k}$. Let
$\widetilde{h}_{\gamma}^{(K)}$ be the online bandwidth given explicitly in
(18) of Section 4 and $\\{\eta_{\gamma,l}^{(K)}\\}_{l=1}^{L}$ be the dynamic
candidate bandwidths at block $K$ satisfying
$\widetilde{h}_{\gamma}^{(K)}=\eta_{\gamma,1}^{(K)}>\ldots\eta_{\gamma,L}^{(K)}$.
Then by setting $\omega_{\gamma,k}^{(K)}=s_{k,2}/S_{K,2}$ and substituting
$\eta_{\mu,l}^{(K)}$ with $\eta_{\gamma,l}^{(K)}$ in (5), (8) and (7), we
obtain the dynamic centroids $\\{\phi_{\gamma,l}^{(k)}\\}_{l=1}^{L}$ and
pseudo-sufficient statistics
$\\{\widetilde{P}_{\gamma,l}^{(K)},\widetilde{q}_{\gamma,l}^{(K)}\\}_{l=1}^{L}$
which are of form
$\\{\sum_{k=1}^{K}P_{k}(s,t;\cdot),\sum_{k=1}^{K}q_{k}(s,t;\cdot)\\}$, where
$P_{k}(s,t;\cdot),q_{k}(s,t;\cdot)$ are defined as in (4). Let
$e_{3}=(1,0,0)^{\mathrm{T}}$, and the online covariance estimate is given by

$\widetilde{\gamma}^{(K)}(s,t)=e_{3}^{\mathrm{T}}\left\\{\widetilde{P}_{\gamma,1}^{(K)}(s,t)\right\\}^{-1}\widetilde{q}_{\gamma,1}^{(K)}(s,t).$
(12)

To conclude, we state the following algorithm for the mean estimation. We
mention that the proposed method can be applied to the general $d$-dimensional
local linear regression as well as higher-order local polynomials. Please see
S.3 of Supplementary Material for more detail.

Algorithm 1 Online Estimation for Mean Function

1: Initialize: $K\leftarrow 0$, $S_{K,1}\leftarrow 0$,
$\phi_{\mu,l}^{(K)}\leftarrow 0$, $\widetilde{P}_{\mu,l}^{(K)}(t)\leftarrow
0$, $\widetilde{q}_{\mu,l}^{(K)}(t)\leftarrow 0$, for $l=1,2,\ldots,L$.

2: While $K<K_{max}$:

3: input: $\\{(T_{Kij},Y_{Kij}):i=1,\ldots,n_{K},j=1,\ldots,m_{Ki}\\}$

4: $K\leftarrow K+1$, $s_{K,1}\leftarrow\sum_{i=1}^{n_{K}}m_{Ki}$,
$S_{K,1}\leftarrow S_{K,1}+s_{K,1}$;

5: compute $\widetilde{h}_{\mu}^{(K)}$ as in (18) based on pilot estimates as
introduced in Appendix B;

6: for $l=1,2,\ldots,L$:

7: compute $\eta_{\mu,l}^{(K)}$ based on (19) in Section 4;

8: compute ${j}_{l}^{(K)}$ based on (8);

9: update $\phi_{\mu,l}^{(K)}$ by (7);

10: update $\widetilde{P}_{\mu,l}^{(K)}(t)$ and
$\widetilde{q}_{\mu,l}^{(K)}(t)$ by (5);

11: end for

12: output:
$\widetilde{\mu}^{(K)}(t)=e_{2}^{\mathrm{T}}\\{\widetilde{P}_{\mu,1}^{(K)}(t)\\}^{-1}\widetilde{q}_{\mu,1}^{(K)}(t)$

13: end while

##  4 Theoretical Analysis

In this section, we illustrate the asymptotic normality of the online mean and
covariance estimates of different sampling schemes. Then we study the
convergence of online bandwidth selection and derive the optimal dynamic
candidate bandwidth sequence, which gives a lower bound for the relative
efficiency versus the batch estimates.

Recall that there are $n_{k}$ subjects in the $k$th block and the $i$th
subject is measured at $m_{ki}$ times for $i=1,\ldots,n_{k}$, and
$s_{K,l},S_{K,l}$ are defined in (6). To describe the phase transition
phenomenon of the mean and covariance estimators, we denote the total number
of subjects up to block $K$ and the individual averages of observations for
$\mu$ and $\gamma$ by,

$N_{K}=\sum_{k=1}^{K}n_{k},\ \overline{m}_{\mu,K}=S_{K,1}/N_{K},\
\overline{m}_{\gamma,K}=S_{K,2}/N_{K}.$ (13)

Write $a\asymp b$ if $a/b\rightarrow C\in(0,\infty)$ as
$a,b\rightarrow\infty$. Recall that $K$ is the total number of blocks, and
$N_{K}\geq K$ is the number of subjects. We stress that the asymptotic results
in this section are for $K\rightarrow\infty$ in spirit of streaming data
accumulation. Recalling that $\omega_{\mu,k}^{(K)}=s_{k,1}/S_{K,1}$ and
$\omega_{\gamma,k}^{(K)}=s_{k,2}/S_{K,2}$, we further define

$\rho_{\mu,j}^{(K)}=\sum_{k=1}^{K}\omega_{\mu,k}^{(K)}\left(\widetilde{\eta}_{\mu,k}^{(K)}\right)^{j},\
\rho_{\gamma,j}^{(K)}=\sum_{k=1}^{K}\omega_{\gamma,k}^{(K)}\left(\widetilde{\eta}_{\gamma,k}^{(K)}\right)^{j},\
j\in\mathbb{Z}.$ (14)

######  Theorem 1.

Under (A.1)–(A.4), (A.6) and (A.8)–(A.9) in Appendix A, let $f$ be the density
of $T$ and $\rho_{\mu,j}^{(K)}$ be defined as in (14). Suppose that
$\widetilde{h}_{\mu}^{(K)}\asymp S_{K,1}^{-1/5}$,
$\widetilde{h}_{\mu}^{(K)}=\eta_{\mu,1}^{(K)}>\cdots>\eta_{\mu,L}^{(K)}$, and
$\lim\sup_{K}\overline{m}_{\gamma,K}/\overline{m}_{\mu,K}^{2}<\infty$, where
$\overline{m}_{\mu,K},\overline{m}_{\gamma,K}$ are defined as in (13). Denote
$R(W)=\int W(x)^{2}dx$, $\alpha(W)=\int x^{2}W(x)dx$ and recall that
$\sigma^{2}$ is the noise variance. For a fixed interior point $t\in(0,1)$, as
$K\rightarrow\infty$, the estimate $\widetilde{\mu}^{(K)}(t)$ as in (9)
satisfies

$\displaystyle\left\\{S_{K,1}\big{/}\rho_{\mu,-1}^{(K)}\right\\}^{\frac{1}{2}}\left\\{\widetilde{\mu}^{(K)}(t)-\mu(t)-\frac{1}{2}\alpha(W)\mu^{\prime\prime}(t)\rho_{\mu,2}^{(K)}+o_{p}\left(\rho_{\mu,2}^{(K)}\right)\right\\}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N\big{(}0,\Gamma_{\mu}(t)\big{)},$

where

$\Gamma_{\mu}(t)=R(W)\frac{\gamma(t,t)+\sigma^{2}}{f(t)}+\gamma(t,t)\frac{S_{K,2}}{S_{K,1}\rho_{\mu,-1}^{(K)}},$

which gives the following phase transition,

  * (1)

when $\overline{m}_{\mu,K}/N_{K}^{1/4}\rightarrow 0$,

$\displaystyle\Gamma_{\mu}(t)=R(W)\frac{\gamma(t,t)+\sigma^{2}}{f(t)};$

  * (2)

when $\overline{m}_{\mu,K}/N_{K}^{1/4}\rightarrow C$, the variance becomes

$\displaystyle\Gamma_{\mu}(t)=R(W)\frac{\gamma(t,t)+\sigma^{2}}{f(t)}+C_{1}\gamma(t,t),$

where $C_{1}=\lim_{K}S_{K,2}/\\{S_{K,1}\rho_{\mu,-1}^{(K)}\\}$, and
$C,C_{1}\in(0,\infty)$;

  * (3)

when $\overline{m}_{\mu,K}/N_{K}^{1/4}\rightarrow\infty$, the bias vanishes
and the distribution is simplified to

$\displaystyle\left\\{N_{K}S_{K,1}^{2}/S_{K,2}\right\\}^{\frac{1}{2}}\Big{\\{}\widetilde{\mu}^{(K)}(t)-\mu(t)\Big{\\}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N\big{(}0,\gamma(t,t)\big{)}.$

Theorem 1 gives a systematic partition of functional data into three
categories according to the numbers of repeated measurements on experimental
subjects, yielding different convergence properties, which shall be discussed
more after representing the asymptotic behaviors of covariance estimate. Note
that in case (2) we have $S_{K,1}/S_{K,2}\asymp\overline{m}_{\mu,K}^{-1}\asymp
S_{K,1}^{-1/5}$ and
$|\rho_{\mu,-1}^{(K)}-1/\widetilde{h}_{\mu}^{(K)}|=|\sum_{k=1}^{K}\omega_{\mu,k}^{(K)}(1/\widetilde{\eta}_{\mu,k}^{(K)})-1/\widetilde{h}_{\mu}^{(K)}|\leq\sum_{k=1}^{K}\omega_{\mu,k}^{(K)}|1/\widetilde{h}_{\mu}^{(k)}-1/\widetilde{h}_{\mu}^{(K)}|$,
then $S_{K,2}/\\{S_{K,1}\rho_{\mu,-1}^{(K)}\\}\rightarrow C_{1}\in(0,\infty)$
holds when $\widetilde{h}_{\mu}^{(K)}\asymp S_{K,1}^{-1/5}$. Though we assume
that $T_{kij}$ are i.i.d. across and within subjects, it is noted that when
$T_{kij}$ may be dependent within subjects, such as scheduled visit times with
random missing and/or fluctuation, the proposed online estimator also attains
the same convergence result as the batch estimator by slightly modifying the
proofs.

Recall that Zhang and Wang (2016) has taken the mean function $\mu$ as known
in the covariance estimation, so $\widehat{\mu}^{(K)}$ is not entangled in the
asymptotic distribution of $\widehat{\gamma}^{(K)}$. In the batch setting, we
use the batch estimator $\widehat{\mu}^{(K)}$ instead and establish the
asymptotic normality of the covariance estimation which is proved to be the
same as using $\mu$. To see this, denote
$e_{kij}^{(K)}=\mu(T_{kij})-\widehat{\mu}^{(K)}(T_{kij})$ and
$\zeta_{kij}=\Phi_{ki}(T_{kij})+\varepsilon_{kij}$, then

$\widehat{C}_{ki,j_{1},j_{2}}=\left(e_{kij_{1}}^{(K)}+\zeta_{kij_{1}}\right)\left(e_{kij_{2}}^{(K)}+\zeta_{kij_{2}}\right),\quad
C_{ki,j_{1},j_{2}}=\zeta_{kij_{1}}\zeta_{kij_{2}}.$

On one hand, $\mathop{\rm Var}(\widehat{C}_{ki,j_{1},j_{2}})$ and $\mathop{\rm
Var}(C_{ki,j_{1},j_{2}})$ are of the same order, both dominated by
$\zeta_{kij_{1}}\zeta_{kij_{2}}$. On the other, note that
$\widehat{\mu}^{(K)}$ is the weighted average of $\\{Y_{kij}:1\leq k\leq
K,1\leq i\leq n_{k},1\leq j\leq m_{ki}\\}$, and $e_{ki,j_{1}}^{(K)}$ is weakly
correlated to $\zeta_{kij_{2}}$, which makes the bias induced by
$e_{ki,j_{1}}^{(K)}\zeta_{kij_{2}}$ negligible. The bias of
$e_{ki,j_{1}}^{(K)}e_{ki,j_{2}}^{(K)}\asymp(\widehat{h}_{\mu}^{(K)})^{4}$ is
also negligible compared to $(\widehat{h}_{\gamma}^{(K)})^{2}$. This argument
also applies to the online estimation, and indicates that the bandwidth
selection for mean and covariance estimators can be decoupled, which greatly
facilitates the derivation. For space economy, we present the result for the
online covariance estimation (using $\widetilde{\mu}^{(K)}$), and defer the
batch version (using $\widehat{\mu}^{(K)}$) to Lemma 1 in Appendix A, while
the detailed proofs are offered in the Supplementary Material. Denote

$\displaystyle V_{1}(s,t)$ $\displaystyle=\mathop{\rm
Var}\big{\\{}\big{(}Y_{1}-\mu(T_{1})\big{)}\big{(}Y_{2}-\mu(T_{2})\big{)}|T_{1}=s,T_{2}=t\big{\\}},$
$\displaystyle V_{2}(s,t)$ $\displaystyle=\mathop{\rm
Cov}\big{\\{}\big{(}Y_{1}-\mu(T_{1})\big{)}\big{(}Y_{2}-\mu(T_{2})\big{)},\big{(}Y_{1}-\mu(T_{1})\big{)}\big{(}Y_{3}-\mu(T_{3})\big{)}|T_{1}=s,T_{2}=t,T_{3}=t\big{\\}},$
(15) $\displaystyle V_{3}(s,t)$ $\displaystyle=\mathop{\rm
Cov}\big{\\{}\big{(}Y_{1}-\mu(T_{1})\big{)}\big{(}Y_{2}-\mu(T_{2})\big{)},\big{(}Y_{3}-\mu(T_{3})\big{)}\big{(}Y_{4}-\mu(T_{4})\big{)}|T_{1}=s,T_{2}=t,T_{3}=s,T_{4}=t\big{\\}}.$

######  Theorem 2.

Under (A.1)–(A.4) and (A.7)–(A.9) in Appendix A, let $\rho_{\gamma,j}^{(K)}$
be defined as in (14). Suppose that $\widetilde{h}_{\gamma}^{(K)}\asymp
S_{K,2}^{-1/6}$,
$\widetilde{h}_{\gamma}^{(K)}=\eta_{\gamma,1}^{(K)}>\cdots>\eta_{\gamma,L}^{(K)}$,
and $f,R(W),\alpha(W),\sigma^{2}$ are the same as in Theorem 1. For a fixed
interior point $(s,t)\in(0,1)^{2}$, as $K\rightarrow\infty$, the covariance
estimator $\widetilde{\gamma}^{(K)}$ in (12) satisfies

$\displaystyle\left\\{S_{K,2}\big{/}\rho_{\gamma,-2}^{(K)}\right\\}^{\frac{1}{2}}\left\\{\widetilde{\gamma}^{(K)}(s,t)-\gamma(s,t)-\frac{1}{2}\alpha(W)\left(\frac{\partial^{2}\gamma}{\partial
s^{2}}+\frac{\partial^{2}\gamma}{\partial
t^{2}}\right)\rho_{\gamma,2}^{(K)}+o_{p}\left(\rho_{\gamma,2}^{(K)}\right)\right\\}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N\big{(}0,\Gamma_{\gamma}(s,t)\big{)},$

where

$\displaystyle\Gamma_{\gamma}(s,t)=$
$\displaystyle\\{1+I(s=t)\\}\left\\{R(W)^{2}\frac{V_{1}(s,t)}{f(s)f(t)}\right.$
$\displaystyle\left.+\frac{1}{S_{K,2}\rho_{\gamma,-2}^{(K)}}\sum_{k=1}^{K}\frac{s_{k,3}}{\widetilde{\eta}_{\gamma,k}^{(K)}}R(W)\frac{f(s)V_{2}(t,s)+f(t)V_{2}(s,t)}{f(s)f(t)}\right\\}+\frac{S_{K,4}}{S_{K,2}\rho_{\gamma,-2}^{(K)}}V_{3}(s,t),$

and the following statements hold,

  * (1)

when $\overline{m}_{\gamma,K}/N_{K}^{1/2}\rightarrow 0$,

$\displaystyle\Gamma_{\gamma}(s,t)=\\{1+I(s=t)\\}\frac{R(W)^{2}V_{1}(s,t)}{f(s)f(t)};$

  * (2)

when $\overline{m}_{\gamma,K}/N_{K}^{1/2}\rightarrow C$,

$\displaystyle\Gamma_{\gamma}(s,t)=$
$\displaystyle\\{1+I(s=t)\\}\left\\{\frac{R(W)^{2}V_{1}(s,t)}{f(s)f(t)}\right.$
$\displaystyle\left.+\frac{R(W)C_{1}}{C_{0}^{1/2}}\frac{f(s)V_{2}(t,s)+f(t)V_{2}(s,t)}{f(s)f(t)}\right\\}+C_{1}^{2}V_{3}(s,t),$

where $C_{0}=\lim_{K}S_{K,2}S_{K,4}/S_{K,3}^{2}$,
$C_{1}^{2}=\lim_{K}S_{K,4}/\\{S_{K,2}\rho_{\gamma,-2}^{(K)}\\}$, and
$C,C_{0},C_{1}\in(0,\infty)$;

  * (3)

when $\overline{m}_{\gamma,K}/N_{K}^{1/2}\rightarrow\infty$, the bias vanishes
and the distribution is simplified to

$\displaystyle\left\\{N_{K}S_{K,2}^{2}\big{/}S_{K,4}\right\\}^{\frac{1}{2}}\Big{\\{}\widetilde{\gamma}^{(K)}(s,t)-\gamma(s,t)\Big{\\}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N\big{(}0,V_{3}(s,t)\big{)}.$

Similar to case (2) of Theorem 1,
$C_{1}^{2}=\lim_{K}S_{K,4}/\\{S_{K,2}\rho_{\gamma,-2}^{(K)}\\}\in(0,\infty)$
holds as long as $\widetilde{h}_{\gamma}^{(K)}\asymp S_{K,2}^{-1/6}$. Note
that the partition in Theorem 2 is identical to that in Theorem 1. Case (1) in
these two theorems is the so-called sparse design, where mean and covariance
estimates attain nonparametric convergence rate that is slower than
root-$N_{K}$. Case (2) and (3) are dense data such that the parametric
convergence rate $N_{K}^{1/2}$ can be achieved. We refer to case (2) as
moderately dense and case (3) ultra dense for that, under the latter scheme,
both the bias term and the variance discontinuity vanish. Such phase
transition phenomenon for online estimation is accordant with the classical
batch result in Zhang and Wang (2016). Using the same techniques, one can
further prove that, in terms of $L_{2}$ and uniform convergences, the online
estimates can also attain the rates as the batch results in Zhang and Wang
(2016).

To derive the optimal bandwidths, we introduce the integrated mean squared
error. For a function $m(x)$ and its estimate $\widehat{m}(x)$, let $f_{X}$ be
the density function of $X$ supported on $I$, then $IMSE(\widehat{m})$ is a
global loss criterion defined as

$IMSE(\widehat{m})=\int_{I}\left[{\mathrm{E}}\\{\widehat{m}(x)-m(x)\\}^{2}\right]f_{X}(x)dx.$

Note that ${\mathrm{E}}\\{\widehat{m}(x)-m(x)\\}^{2}=\mathop{\rm
Var}\widehat{m}(x)+[{\mathrm{E}}\\{\widehat{m}(x)\\}-m(x)]^{2}$, we have the
following results for the batch and online estimates.

######  Corollary 1.

Under (A.1)–(A.4) and (A.6)–(A.9) in Appendix A, as $K\rightarrow\infty$,

$\displaystyle IMSE\left(\widehat{\mu}^{(K)}\right)$
$\displaystyle=\frac{1}{4}\theta_{\mu}\alpha^{2}(W)\big{\\{}\widehat{h}_{\mu}^{(K)}\big{\\}}^{4}+\frac{\nu_{\mu}}{S_{K,1}\widehat{h}_{\mu}^{(K)}}+o_{p}\left(\big{\\{}\widehat{h}_{\mu}^{(K)}\big{\\}}^{4}+S_{K,1}^{-1}\big{\\{}\widehat{h}_{\mu}^{(K)}\big{\\}}^{-1}\right),$
$\displaystyle IMSE\left(\widehat{\gamma}^{(K)}\right)$
$\displaystyle=\frac{1}{4}\theta_{\gamma}\alpha^{2}(W)\big{\\{}\widehat{h}_{\gamma}^{(K)}\big{\\}}^{4}+\frac{\nu_{\gamma}}{S_{K,2}\big{\\{}\widehat{h}_{\gamma}^{(K)}\big{\\}}^{2}}+o_{p}\left(\big{\\{}\widehat{h}_{\gamma}^{(K)}\big{\\}}^{4}+S_{K,2}^{-1}\big{\\{}\widehat{h}_{\gamma}^{(K)}\big{\\}}^{-2}\right),$

where

$\displaystyle\theta_{\mu}=\int\mu^{\prime\prime}(t)f(t)dt,\
\theta_{\gamma}=\int\int(\partial^{2}\gamma/\partial
s^{2}+\partial^{2}\gamma/\partial t^{2})f(s)f(t)dsdt,$
$\displaystyle\nu_{\mu}=\int\Gamma_{\mu}(t)f(t)dt,\
\nu_{\gamma}=\int\int\Gamma_{\gamma}(s,t)f(s)f(t)dsdt.$ (16)

######  Corollary 2.

Suppose that the conditions in Theorem 1–2 hold, then as $K\rightarrow\infty$,

$\displaystyle IMSE\left(\widetilde{\mu}^{(K)}\right)$
$\displaystyle=\frac{1}{4}\theta_{\mu}\alpha^{2}(W)\big{\\{}\rho_{\mu,2}^{(K)}\big{\\}}^{2}+\frac{\nu_{\mu}}{S_{K,1}}\rho_{\mu,-1}^{(K)}+o_{p}\left(\left\\{\rho_{\mu,2}^{(K)}\right\\}^{2}+S_{K,1}^{-1}\rho_{\mu,-1}^{(K)}\right),$
$\displaystyle IMSE\left(\widetilde{\gamma}^{(K)}\right)$
$\displaystyle=\frac{1}{4}\theta_{\gamma}\alpha^{2}(W)\big{\\{}\rho_{\gamma,2}^{(K)}\big{\\}}^{2}+\frac{\nu_{\gamma}}{S_{K,2}}\rho_{\gamma,-2}^{(K)}+o_{p}\left(\left\\{\rho_{\gamma,2}^{(K)}\right\\}^{2}+S_{K,2}^{-1}\rho_{\gamma,-2}^{(K)}\right),$

where $\theta_{\mu},\theta_{\gamma},\nu_{\mu},\nu_{\gamma}$ are defined as in
Corollary 1.

From Corollary 1, the optimal bandwidths to minimize
$IMSE\left(\widehat{\mu}^{(K)}\right)$ and
$IMSE\left(\widehat{\gamma}^{(K)}\right)$ are, respectively,

$h_{\mu,*}^{(K)}=\left(\frac{\nu_{\mu}}{\alpha^{2}(W)\theta_{\mu}}\right)^{\frac{1}{5}}S_{K,1}^{-\frac{1}{5}},\quad
h_{\gamma,*}^{(K)}=\left(\frac{\nu_{\gamma}}{\alpha^{2}(W)\theta_{\gamma}}\right)^{\frac{1}{6}}S_{K,2}^{-\frac{1}{6}}.$
(17)

We suggest to use pilot estimates adopting the online method to approximate
the unknown integrals $\theta_{\mu},\theta_{\gamma}$ and the variance terms
$\nu_{\mu},\nu_{\gamma}$, see details in Appendix B. The online estimates of
optimal bandwidths are

$\widetilde{h}_{\mu}^{(K)}=\left(\frac{\widetilde{\nu}_{\mu}^{(K)}}{\alpha^{2}(W)\widetilde{\theta}_{\mu}^{(K)}}\right)^{\frac{1}{5}}S_{K,1}^{-\frac{1}{5}},\
\widetilde{h}_{\gamma}^{(K)}=\left(\frac{\widetilde{\nu}_{\gamma}^{(K)}}{\alpha^{2}(W)\widetilde{\theta}_{\gamma}^{(K)}}\right)^{\frac{1}{6}}S_{K,2}^{-\frac{1}{6}},$
(18)

where $\widetilde{\theta}_{\mu}^{(K)},\widetilde{\theta}_{\gamma}^{(K)}$ and
$\widetilde{\nu}_{\mu}^{(K)},\widetilde{\nu}_{\gamma}^{(K)}$ are given
explicitly in Appendix B. Let $\widehat{h}_{\mu}^{(K)}$ be the bandwidth
selected by replacing the online pilot estimates in (18) with the batch ones.
When the bandwidths for pilot estimates are of appropriate orders, the online
bandwidths (18) can attain the same convergence rates as the batch
competitors, as stated in the following theorem.

######  Theorem 3.

Suppose that the conditions in Theorem 1–2 hold. Let
$\widetilde{h}_{\mu}^{(K)}$ and $\widetilde{h}_{\gamma}^{(K)}$ be the proposed
online bandwidth estimates (18) where
$\widetilde{\theta}_{\mu}^{(K)},\widetilde{\theta}_{\gamma}^{(K)}$ and
$\widetilde{\nu}_{\mu}^{(K)},\widetilde{\nu}_{\gamma}^{(K)}$ are based on
online local polynomials with bandwidths and candidates given explicitly in
(B.1) of Appendix B. Then the online bandwidth selection satisfies that as
$K\rightarrow\infty$,

$\frac{\widetilde{h}_{\mu}^{(K)}-h_{\mu,*}^{(K)}}{h_{\mu,*}^{(K)}}=O_{p}\left(S_{K,1}^{-\frac{2}{7}}\right),\quad\frac{\widetilde{h}_{\gamma}^{(K)}-h_{\gamma,*}^{(K)}}{h_{\gamma,*}^{(K)}}=O_{p}\left(S_{K,2}^{-\frac{1}{4}}\right),$

where $h_{\mu,*}^{(K)},h_{\gamma,*}^{(K)}$ are the optimal bandwidths (17).

Define the relative efficiency in terms of integrated mean squared errors by

$eff_{K}(\widetilde{\mu}^{(K)})=IMSE(\widetilde{\mu}^{(K)})/IMSE(\widehat{\mu}^{(K)}),\quad
eff_{K}(\widetilde{\gamma}^{(K)})=IMSE(\widetilde{\gamma}^{(K)})/IMSE(\widehat{\gamma}^{(K)}),$

which measures the performance of the proposed online method compared to the
classical estimate using full data. We now present our key result, the choice
of candidate bandwidth sequence which makes $\widetilde{\eta}_{\mu,k}^{(K)}$
(or $\widetilde{\eta}_{\gamma,k}^{(K)}$) as close to
$\widetilde{h}_{\mu}^{(K)}$ (or $\widetilde{h}_{\gamma}^{(K)}$) as possible
for all $k$ and maximizes the relative efficiency with a lower bound.

######  Theorem 4.

Under (A.1)–(A.9) in Appendix A, when the candidate bandwidths are

$\eta_{\mu,l}^{(K)}=\left(\frac{L-l+1}{L}\right)^{\frac{1}{5}}\widetilde{h}_{\mu}^{(K)},\quad\eta_{\gamma,l}^{(K)}=\left(\frac{L-l+1}{L}\right)^{\frac{1}{6}}\widetilde{h}_{\gamma}^{(K)},\quad
l=1,\ldots,L,$ (19)

then as $K\rightarrow\infty$, the corresponding online estimates
$\widetilde{\mu}^{(K)}$ and $\widetilde{\gamma}^{(K)}$ attain the optimal
relative efficiency with the following lower bound,

$\left\\{1+\left(\frac{2d}{d+4}c_{1}+\frac{4}{d+4}c_{2}\right)\frac{1}{L}+\frac{d}{d+4}c_{2}^{2}\frac{1}{L^{2}}\right\\}^{-1}+O_{p}\left(S_{K,d}^{-\frac{2}{d+6}}\right),$

where $d=1$ for $\widetilde{\mu}^{(K)}$ and $d=2$ for
$\widetilde{\gamma}^{(K)}$. Specifically,

$\displaystyle
eff(\widetilde{\mu}^{(K)})\geq(1+0.1831L^{-1}+0.0032L^{-2})^{-1}+O_{p}\left(S_{K,1}^{-2/7}\right),$
$\displaystyle
eff(\widetilde{\gamma}^{(K)})\geq(1+0.2422L^{-1}+0.0190L^{-2})^{-1}+O_{p}\left(S_{K,2}^{-1/4}\right).$
(20)

From the illustration in the left panel of Figure 4, the relative efficiencies
of the proposed online mean and covariance estimates increase rapidly to
exceed 95% when $L\geq 5$. This is desirable to attain high efficiency with a
small/moderate $L$. Note that the computational cost in terms of time and
memory is proportional to $L$, this lower bound helps make an informed trade-
off between statistical and computational efficiency, which makes the proposed
method practically useful. We also define the relative error,

$\delta(\widetilde{\mu}^{(K)})=1-eff(\widetilde{\mu}^{(K)}),\quad\delta(\widetilde{\gamma}^{(K)})=1-eff(\widetilde{\gamma}^{(K)}),$

whose upper bound can be derived from (4). We plot $L$ versus the upper bound
of $\log(\delta)$ in the right panel of Figure 4 as a suggestion for the
selection of $L$ in practice.

Figure 4:  The left panel shows lower bounds for the relative efficiencies of
the proposed online mean (solid lines) and covariance (dashed lines) estimates
versus different lengths $L$ of candidate bandwidth sequences. The right panel
plots the selection of $L$ under different upper bounds of $\log(\delta)$ for
the proposed online mean (solid lines) and covariance (dashed lines)
estimates.

We conclude this section by mentioning that the asymptotic distributions of
$\widetilde{\mu}^{(K)}$ and $\widetilde{\gamma}^{(K)}$ have the consistent
form with the general $d$-dimensional local linear regression. Hence Theorem 3
and 4 hold for general local linear regression with $d$-dimensional
covariates. See details in S.4 of Supplementary Material.

##  5 Simulation

We conduct simulation to illustrate the performance of the proposed online
method and verify the theoretical findings in Section 4. Let noises
${\varepsilon}_{kij}$ be i.i.d. from $N(0,0.5^{2})$ and $T_{kij}$ follow a
uniform distribution on $[0,1]$. Data are generated by
$Y_{kij}=X_{ki}(T_{kij})+{\varepsilon}_{kij}$, where
$X_{ki}(T_{kij})=\mu(T_{kij})+\Phi_{ki}(T_{kij})$ with the mean function
$\mu(t)=2\sin(2\pi t)$ and the stochastic part
$\Phi(t)=\sum_{i=1}^{10}\xi_{i}\phi_{i}(t)$, where $\phi_{1}(t)=1$,
$\phi_{i}(t)=\sqrt{2}\cos\\{(i-1)\pi t\\}$, $\xi_{i}$ are independently
sampled from $N(0,\lambda_{i})$ with $\lambda_{i}=0.4\times i^{-2},\
i=1,\ldots,10$. Recall that there are $n_{k}$ subjects in the $k$th block,
among which the $j$th one has $m_{ki}$ observations. For sparse data, we let
$n_{k}$ follow the normal distribution $N(20,9)$ and $m_{ki}$ follow $N(6,4)$;
for dense data, we let $n_{k}=3$ and $m_{ki}$ follow $N(20,4)$; which are all
rounded off to the nearest integers. The experiment is repeated 100 times,
each with $K_{max}=1000$ blocks.

We derive in (B.1) of Appendix B that bandwidths for pilot estimates of
$\theta_{\mu},\theta_{\gamma}$ and $\nu_{\mu},\nu_{\gamma}$ in (1) shall be

$h_{\theta_{\mu}}^{(K)}=GS_{K,1}^{-1/7},\
h_{\nu_{\mu}}^{(K)}=RS_{K,1}^{-1/5},\
h_{\theta_{\gamma}}^{(K)}=GS_{K,2}^{-1/8},\
h_{\nu_{\gamma}}^{(K)}=RS_{K,2}^{-1/6}.$

After experimenting an extensive range of $G$ and $R$, we find that the
convergence rate of the online bandwidth is not sensitive to the values of
these two quantities. Values in $[0.5^{1/d},1]$ are in general adequate, thus
we set $G=R=0.5^{1/d}$ in the sequel, where $d=1$ for estimating $\mu$ and
$d=2$ for estimating $\gamma$. Let $J$ be the length of candidate bandwidth
sequences for the pilot estimates of $\theta_{\mu},\theta_{\gamma}$ and
$\nu_{\mu},\nu_{\gamma}$, and recall that $L$ is the length of candidate
bandwidth sequences for estimating $\mu$ and $\gamma$, we set $J=L$ when
estimating $\mu$ and $J=3$ when estimating $\gamma$ to ease computation. We
also stop updating $\widetilde{\theta}_{\gamma}^{(K)}$ and
$\widetilde{\nu}_{\gamma}^{(K)}$ after $K^{\prime}=200$ to further reduce
computation without influencing the convergence rate as
$K^{\prime}/K_{max}=O_{p}(1)$. Note that the classical batch method is
computational expensive when sample size tends large, we implement it at every
40 blocks, i.e. $K=40,80,120,\ldots$. For mean and covariance estimation of
sparse and dense data, we examine the following measures.

1\. Relative efficiency. Figure 5 shows the empirical relative efficiencies of
the mean and covariance estimates that increase with $L$ and are stably higher
than the theoretical lower bounds in Theorem 4 when $K$ tends large. Recall
that Kong and Xia (2019) studied the $d$-dimensional online local linear
regression of independent noises, which is not applicable to functional data.
Even under the assumption of independent errors, their theoretical lower
bounds of relative efficiency, which are 0.943 when $d=1$ and 0.924 when
$d=2$, are still smaller than ours. An explanation is that their stored
statistics are fixed, and if the bandwidths of previous blocks deviate from
the optimal values, there is no data-driven adjustment in the sequential
estimates. By contrast, our stored statistics are updated dynamically by
selecting different bandwidths from the candidate sequences and hence produce
more efficient estimates.

Figure 5:  The empirical relative efficiencies of the proposed online mean and
covariance estimates for functional data based on sparse design (solid) and
dense design (dashed) versus the theoretical lower bound (dotted) with
candidate bandwidth sequence length $L=3,5,10$ and $20$.

2\. Bandwidth selection. The convergence of bandwidths in Theorem 3 is also
examined in Figure 6, where both batch and online selections converge to the
theoretical optimal bandwidth along with data collection. Moreover, we depict
the dynamically updated pseudo-bandwidth sequences
$\\{\widetilde{\eta}_{\mu,k}^{(K)}\\}_{k=1}^{K}$ and
$\\{\widetilde{\eta}_{\gamma,k}^{(K)}\\}_{k=1}^{K}$ that are used to produce
the estimates at time $K=200,500,1000$, similar to that shown in Figure 3.
This provides empirical support for the fact that the our method is indeed
able to adjust the sub-statistics in previous blocks in spite of no access to
those data.

Figure 6:  The four rows show, respectively, the bandwidth selection for mean
and covariance estimation of sparse and dense data. The first column shows the
Monto Carlo averages of the bandwidths selected by the online method (solid)
$\widetilde{h}_{\mu}^{(K)},\widetilde{h}_{\gamma}^{(K)}$ and the batch method
$\widehat{h}_{\mu}^{(K)},\widehat{h}_{\gamma}^{(K)}$ (dashed) based on 100
runs for $K=1,2,\ldots,1000$, and the other three columns depict the
dynamically updated bandwidths
$\\{\widetilde{\eta}_{\mu,k}^{(K)}\\}_{k=1}^{K},\\{\widetilde{\eta}_{\gamma,k}^{(K)}\\}_{k=1}^{K}$
(thick dots) at time $K=200,500,1000$, respectively, along with the dynamic
candidate sequences (light dots) and optimal bandwidth estimates (connected by
solid line) at each $k$.

3\. Computational time. We compare the computing times of the classical batch
and the proposed online methods using the Unix server of 2.10GHz CPU and 188G
memory with 176 logic cores. It is noted from Figure 7 that the computing time
of the classical batch method grows approximately linearly with data blocks,
while the online method spends nearly constant time given similar block sizes.
Thus we graphed only the first 400 blocks for visualization when compared with
the online method using $L=3,5,10,20$ that achieves substantial computational
saving proportional to $L$.

Figure 7:  The comparison of computing times between the proposed online
method and the batch method when estimating mean and covariance functions
based on sparse and dense functional data. The dashed lines correspond to the
computing times of the batch method and the solid lines to the online method,
from the lowest to the highest in each panel representing $L=3,5,10,20$,
respectively.

We close this section by suggesting a reasonable range of $L\in[5,20]$ based
on our empirical and theoretical findings. The choice shall be made depending
on whether accuracy or computation is of main concern.

##  6 Real Data Examples

In this section, we present two real data examples to illustrate the
usefulness of the proposed online method.

###  6.1 Airline delay example

The airline dataset consists of flight arrival and departure details for all
commercial airports in the USA, from January 1989 to December 2000
(https://community.amstat.org/jointscsg-section/dataexpo/dataexpo2009). This
involves 304 airports whose daily numbers of flights range from 1 to 1269. We
are interested in the pattern of delay time (unit: minute) during peak hours,
i.e., from 6:00 to 23:00. Airports with flights less than 50 per day are
removed. For each airport, flights departed during off-peak hours are also
removed. Suppose that the departure delay time is independent across days and
we treat the daily airport as a subject with flight departure delay times as
measurements. Let one block represent a day and there are 4383 blocks in
total. We randomly select $n_{k}$ airports from the $k$th day to form the
$k$th data block, where $n_{k}$ are discretely distributed with equal
probability among $\\{6,7,\ldots,15\\}$. To obtain comparable results, the
dense and sparse designs use the same airports for each block, and the dense
measurements contain the sparse ones for each subject. Recall the number of
measurements of the $j$th subject in the $k$th block is denoted by $m_{ki}$.
For dense scheme, we let $m_{ki}$ follow the discrete uniform distribution on
$\\{15,16,\ldots,20\\}$ and randomly pick $m_{ki}$ flights with equal
probability as measurements; for sparse case, choose $m_{ki}$ equal-likely
from 8 to 10 to form the measurements, where $j=1,\ldots,n_{k}$,
$k=1,\ldots,4383$. We also plot the results obtained from all airports and all
flights as baseline for comparison.

We set $L=J=10$ for the estimation of mean function, and $L=10$ and $J=3$ for
covariance estimation, where $L$ is the length of candidate bandwidth sequence
for main estimation, and $J$ is for the pilot bandwidths. As in simulation, we
update the pilot estimates in bandwidth selection for covariance estimation
till $K=200$. Delay aggravates along time, which is probably due to the
cumulative effect of congestion for flight runways. For both designs, the
online and the batch estimates agree well when more data enter the model. We
also conduct functional principal component analysis based on the covariance
estimates. We show the first two eigenfunctions in Figure 9 at time
$K=365,1461,4383$ (the end of year 1989, 1992, 2000), whose aggregated
fractions of variation explained (FVE) exceeds 95%. In both dense and sparse
designs, as data accumulate, the online estimates become closer to the batch
ones as well as the baseline results using all airports and flights. The first
component has similar pattern to the mean function showing more variation in
later hours, while the second represents the influence of different handling
strategies of airports. An explanation may be that those airports with sound
management system and adequate resource matching can alleviate the aggravating
effect of delay to a certain extent, and vice versa. Comparison of computation
times is given in Figure 10, where the online method shows substantial gains
over the batch method.

Figure 8: Mean and covariance estimates for the airline delay dataset at
$K=365,1461,4383$ (the end of year 1989,1992,2000), where the dashed lines and
light-colored surfaces represent the batch results, and the solid lines and
dark-colored surfaces represent the proposed online estimates, with color blue
representing the dense case and red representing the sparse case. The
estimates based on all airports and flights are also plotted in gray as the
baseline.  Figure 9: Estimates of the first two eigenfunctions
$\phi_{1},\phi_{2}$ for the airline delay dataset at $K=365,1461,4383$ (the
end of year 1989, 1992, 2000), where the dashed lines correspond to the batch
results and the solid to the proposed online method with dense design plotted
in blue and sparse in red. The first eigenfunction $\phi_{1}$ has similar
pattern to the mean function, and occupies 88% FVE. The estimates based on all
airports and flights are also plotted in gray as the baseline. Figure 10:  The
comparison of computing times for the airline delay dataset between our method
and the batch method of the first 1821 blocks/five years for mean and
covariance estimation. The dashed lines correspond to the time complexity of
the batch method and the solid lines to our method.

###  6.2 Online news example

The second example is the online news data, which contains 93239 news items
published during 2015/11/08 to 2016/07/07 and their hourly social feedback
(number of clicks) of 48 hours after publishing on LinkedIn, Facebook and
Google, which can be downloaded from
http://archive.ics.uci.edu/ml/datasets/News+Popularity+in+Multiple+Social+Media+Platforms?tdsourcetag=s_pctim_aiomsg.
Denote the publish time for the $j$th news as $t_{0,j}$ and measurement time
as $t_{ij}=t_{0,j}+i\Delta t$, where $i=1,\ldots,48$, $t_{0,j}\in(0,24]$ and
$\Delta t=1$ (unit: hour). The release time $t_{0,j}$ can be viewed as random,
and hence the observations are irregular by space. We refer to the news
published during 0:00 to 6:00 as morning news and aim to study their number of
clicks per hour over the domain $[6,30]$. i.e., 24 hours after 6:00. We add up
the number of clicks from the four platforms for the morning news and further
remove the inactive news whose total number of clicks is less than 6. The
remaining 5689 news are divided into 1138 blocks. We set the blockwise number
of subjects as $n_{k}=5$ for $k\leq 1137$ and $n_{k}=4$ for $k=1138$. For each
subject, let the dense measurements contain the sparse ones. For dense scheme,
we set $m_{ki}=24$ and for sparse scheme $m_{ki}=5$ with measurements randomly
chosen from the dense set with equal probability, $j=1,\ldots,n_{k}$,
$k=1,\ldots,1138$. Let $Y_{kij}^{0}$ denote the original hourly click number,
and we make the following transformation: $Y_{kij}=\log(Y_{kij}^{0}+1)$.

Parameter $L$ and $J$ are set the same as in the airline delay data and the
pilot estimates in bandwidth for covariance estimation update till $K=200$.
Estimates for the mean and covariance function, and eigenfunctions are shown
in Figure 11-12, where difference between online estimates and batch ones
narrows down when data accumulate. Mean function estimates in Figure 11
reveals that the hourly clicks attain peaks at 14:00 and 00:00, i.e., after
the lunch break and before night sleeping. Here we plot the first four
eigenfunction estimates whose aggregated FVE is about 95%. The first
eigenfunction indicates again that the attention a news received is influenced
by people’s daily activity pattern. The second represents sustained growth or
decrease of the news consumption, which is determined by the news quality, as
people like to share valuable news, while articles that attract hits by
headlines are gradually ignored. Finally, one can conclude from Figure 13 that
our online method saves substantial computational time compared to the batch
method.

Figure 11: Mean and covariance estimates for the online news dataset when
$K=100,200,1138$, where the dashed lines and light-colored surfaces correspond
to the batch results, and the solid lines and light-colored surfaces to the
online estimates, with color blue representing the dense case and red
representing the sparse case. The first eigenfunction $\phi_{1}$ has similar
pattern to the mean function, and occupies nearly 66% FVE. Figure 12:
Estimates of the first four eigenfunctions $\phi_{1},\ldots,\phi_{4}$ for the
online news dataset using the online (solid) and batch (dashed) methods at
time $K=100,200,1138$ with dense design plotted in blue and sparse in red.
Figure 13: The comparison of computing times for the online news dataset
between our method and the batch method from $K=1$ to $K=1138$ for mean and
covariance estimation. The dashed lines correspond to the time complexity of
the batch method and the solid lines to our method.

##  7 Concluding Remarks

In this work, we propose a dynamic candidate bandwidth method and apply it to
functional data analysis for the mean and covariance estimation in the online
context. By pre-storing a sequence of statistics and selecting appropriate
candidates, our method dynamically updates statistics for all blocks. To the
best of our knowledge, we are the first to implement dynamic update of
previous statistics with changing bandwidths for streaming data. The proposed
method is both computational and statistical efficient. It requires nearly
constant memory and computing time and has a lower bound for the relative
efficiency in terms of integrated mean squared errors. The computation and
estimation trade-off is achieved by the parameter $L$ that provides useful
guidance in practice.

The proposed method can be adopted to the higher-order local polynomials as
well as other kernel-type estimates. For instance, one can combine the dynamic
candidate bandwidths with profiling algorithm to fit semiparametric partially
linear models, or with backfitting algorithm for additive models. Due to the
interaction between component estimates in the iterative algorithm, the
pseudo-sufficient statistics for each component function would depend on the
estimates of other components and the relative efficiency needs to be further
investigated.

## Supplementary Material

The supplementary material contains the proofs for Theorem 1-4 and Lemma 1.

## Appendix

##  A Assumptions and Auxiliary results

Recall that we observe $Y_{kij}=X_{ki}(T_{kij})+{\varepsilon}_{kij}$, where
$X_{ki}(T_{kij})=\mu(T_{kij})+\Phi_{ki}(T_{kij})$,
$j=1,\ldots,m_{ki},i=1,\ldots,n_{k},k=1,\ldots,K$, with $\mathop{\rm
Cov}\\{X(s),X(t)\\}=\gamma(s,t)$. Denote
$\overline{m}_{j,K}=S_{K,j}/N_{K}\textrm{ for }j\geq 3$. The following
assumptions are imposed for theoretical analysis in Section 4.

  1. (A.1)

The observed time points $\\{T_{kij}\\}$ are i.i.d. copies of a random
variable $T$ defined on [0,1]. The density $f(\cdot)$ of $T$ is bounded away
from 0 and its second derivative is bounded.

  2. (A.2)

$X$ is independent of $T$, and ${\varepsilon}$ is independent of $T$ and
$\Phi$.

  3. (A.3)

The second and fourth derivatives of $\mu$ and $\gamma$ are bounded and
continuous.

  4. (A.4)

The sample size satisfies
$\lim_{K}\overline{m}_{\gamma,K}\cdot\overline{m}_{4,K}/\overline{m}_{3,K}^{2}=C_{0}\in[1,\infty]$,
$\lim\sup_{K}\overline{m}_{\gamma,K}/\overline{m}_{\mu,K}^{2}$,
$\lim\sup_{K}\overline{m}_{\mu,K}\cdot\overline{m}_{3,K}/\overline{m}_{\gamma,K}^{2}$,
$\lim\sup_{K}\overline{m}_{4,K}/\overline{m}_{\gamma,K}^{2}<\infty$, where
$\overline{m}_{l,K},\overline{m}_{\mu,K},\overline{m}_{\gamma,K}$ are defined
in (13).

  5. (A.5)

$\max_{k\leq K}s_{k,j}/S_{K,j}\rightarrow 0$ as $K\rightarrow\infty$ for
$j=1,2$.

  6. (A.6)

$\sup_{t\in[0,1]}{\mathrm{E}}|\Phi(t)|^{3}<\infty$ and
${\mathrm{E}}|{\varepsilon}|^{3}<\infty$.

Assumption (A.1)–(A.4) are general for functional mean and covariance
estimation, see Zhang and Wang (2016). Assumption (A.5) is reasonable to
require the number of measurements in each block is relatively small compared
to the total number of measurements in all $K$ blocks. Assumption (A.6) is a
standard requirement for mean estimation. For covariance estimation, further
assumption on the moments of the random part is needed.

  1. (A.7)

$\sup_{t\in[0,1]}{\mathrm{E}}|\Phi(t)|^{6}<\infty$ and
${\mathrm{E}}|{\varepsilon}|^{6}<\infty$.

We also impose the following assumptions on the kernel function.

  1. (A.8)

$W(\cdot)$ is a symmetric probability density function on $[-1,1]$ and
$\alpha(W),R(W)<\infty$.

  2. (A.9)

$W(\cdot)$ is Lipschitz continuous: There exists $0<L<\infty$ such that

$|W(u)-W(v)|\leq L|u-v|,\ \ {\rm for\ any\ }u,v\in[0,1].$

The lemma below states the asymptotic normality of the batch covariance
estimate based on the batch estimator $\widehat{\mu}^{(K)}$, corresponding to
Theorem 2, which also improves the result of Zhang and Wang (2016) that used
the underlying mean function $\mu$.

######  Lemma 1.

Under (A.1)–(A.4) and (A.7)–(A.9) in Appendix A, suppose that
$\widehat{h}_{\gamma}^{(K)}-h_{\gamma,*}^{(K)}=o_{p}(N_{K}^{-1/6})$, and let
$f$ be the density of $T$. Denote $R(W)=\int W(x)^{2}dx$, $\alpha(W)=\int
x^{2}W(x)dx$ and recall that $\sigma^{2}$ is the noise variance. For a fixed
interior point $(s,t)\in(0,1)^{2}$, as $K\rightarrow\infty$, the covariance
estimator $\widehat{\gamma}^{(K)}$ in (4) satisfies

$\displaystyle\left\\{S_{K,2}\left(\widehat{h}_{\gamma}^{(K)}\right)^{2}\right\\}^{\frac{1}{2}}\left\\{\widehat{\gamma}^{(K)}(s,t)-\gamma(s,t)-\frac{1}{2}\alpha(W)\left(\frac{\partial^{2}\gamma}{\partial
s^{2}}+\frac{\partial^{2}\gamma}{\partial
t^{2}}\right)\left(\widehat{h}_{\gamma}^{(K)}\right)^{2}+o_{p}\left(\left(\widehat{h}_{\gamma}^{(K)}\right)^{2}\right)\right\\}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N\big{(}0,\Gamma_{\gamma}^{\prime}(s,t)\big{)},$

where

$\displaystyle\Gamma_{\gamma}^{\prime}=$
$\displaystyle\\{1+I(s=t)\\}\left\\{S_{K,2}^{-1}R(W)^{2}\frac{V_{1}(s,t)}{f(s)f(t)}\right.$
$\displaystyle\left.+\widehat{h}_{\gamma}^{(K)}S_{K,2}^{-2}S_{K,3}R(W)\frac{f(s)V_{2}(t,s)+f(t)V_{2}(s,t)}{f(s)f(t)}\right\\}+V_{3}(s,t)\left(\widehat{h}_{\gamma}^{(K)}\right)^{-2}S_{K,4}/S_{K,2}^{2},$

The proof for the above lemma is similar to the online case which can be found
in S.2 of the Supplementary Material.

##  B Bandwidth Selection

Recall that the optimal bandwidths for estimating $\mu$ and $\gamma$ are,
respectively,

$\displaystyle
h_{\mu,*}^{(K)}=\left(\frac{\nu_{\mu}}{\alpha^{2}(W)\theta_{\mu}}\right)^{\frac{1}{5}}S_{K,1}^{-\frac{1}{5}},\quad
h_{\gamma,*}^{(K)}=\left(\frac{\nu_{\gamma}}{\alpha^{2}(W)\theta_{\gamma}}\right)^{\frac{1}{6}}S_{K,2}^{-\frac{1}{6}},$

where $\theta_{\mu}=\int\mu^{\prime\prime}(t)f(t)dt$,
$\theta_{\gamma}=\int\int(\partial^{2}\gamma/\partial
s^{2}+\partial^{2}\gamma/\partial t^{2})f(s)f(t)dsdt$,
$\nu_{\mu}=\int\Gamma_{\mu}(t)f(t)dt$ and
$\nu_{\gamma}=\int\int\Gamma_{\gamma}(s,t)f(s)f(t)dsdt$, where
$\Gamma_{\mu},\Gamma_{\gamma}$ are delineated below. We suggest to use pilot
estimates adopting online approach to approximate these quantities. Take
$\theta_{\mu}$ as an example. We first estimate $\mu^{\prime\prime}$ by an
online local cubic smoother with candidate bandwidth sequence
$\\{\eta_{\theta_{\mu},j}^{(k)}\\}_{j=1}^{J}$, $k=1,2,\ldots,K$, and use
$\widetilde{f}^{(K)}=e_{2}^{\mathrm{T}}\widetilde{P}_{\theta_{\mu},1}^{(K)}e_{2}$
to estimate $f$, then integrate
$(\widetilde{\mu}^{\prime\prime})^{(K)}\widetilde{f}_{K}$ to obtain
$\widetilde{\theta}_{\mu}^{(K)}$. And $\theta_{\nu}$ can be estimated by the
same method, with candidate sequence denoted by
$\\{\eta_{\theta_{\gamma},j}^{(k)}\\}_{j=1}^{J}$. When estimating $\nu_{\mu}$
and $\nu_{\gamma}$, the density function can be estimated by
$e_{2}^{\mathrm{T}}\widetilde{P}_{\nu_{\mu},1}^{(K)}e_{2}$ and
$e_{3}^{\mathrm{T}}\widetilde{P}_{\nu_{\gamma},1}^{(K)}e_{3}$, and we need
only estimate $\Gamma_{\mu}$ and $\Gamma_{\gamma}$, which have distinguished
expressions in the cases listed below and require different estimation
procedures. These estimates are based on other local polynomials whose
bandwidth selection is specified after presenting the methodology.

  1. (a)

Mean function estimation for sparse data. Recall the definitions of
$\gamma(t,t)$ and $\sigma^{2}$ in the underlying models are (1), and

$\displaystyle\Gamma_{\mu}=R(W)\frac{\gamma(t,t)+\sigma^{2}}{f(t)}.$

We need only estimate $r(t)=\gamma(t,t)+\sigma^{2}$. We follow the method of
Fan and Yao (1998) to estimate $r(t)$:

     * 1.

Estimate $\mu(t)$ by an online local linear smoother denoted by
$\check{\mu}^{(K)}$ with bandwidth $h_{\nu_{\mu}}^{(K)}=O_{p}(S_{K,1}^{-1/5})$
and the corresponding candidate bandwidth sequence
$\\{\eta_{\nu_{\mu},j}^{(K)}\\}_{j=1}^{J}$, $k=1,2,\ldots,K$;

     * 2.

Compute $\widehat{r}_{Kij}=\\{Y_{Kij}-\check{\mu}^{(K)}(T_{Kij})\\}^{2}$, and
smooth $\widehat{r}_{Kij}$ by an online local linear smoother with the above
$h_{\nu_{\mu}}^{(K)}$ and $\\{\eta_{\nu_{\mu},j}^{(K)}\\}_{j=1}^{J}$ to obtain
$\widetilde{r}^{(K)}(t)$.

  2. (b)

Covariance function estimation for sparse data. Recall the definition of
$V_{1}$ in (4) and

$\displaystyle\Gamma_{\gamma}=\\{1+I(s=t)\\}\frac{R(W)^{2}V_{1}(s,t)}{f(s)f(t)}.$

We only need estimate $V_{1}(s,t)$ that equals
$V_{1}(s,t)={\mathrm{E}}\Phi(s)^{2}\Phi(t)^{2}+\sigma^{2}\\{\gamma(s,s)+\gamma(t,t)\\}+\sigma^{4}-\gamma(s,t)^{2}$,
where $\sigma^{2}$ is the noise variance and $\Phi$ is the stochastic part of
$X$. Estimating $V_{1}(s,t)$ consists of three parts:

     * 1.

Estimate $\gamma(s,t)$ by an online local linear smoother with bandwidth
$h_{\nu_{\gamma}}^{(K)}=O_{p}(S_{K,2}^{-1/6})$ and the corresponding candidate
bandwidth sequence $\\{\eta_{\nu_{\gamma},j}^{(k)}\\}_{j=1}^{J}$,
$k=1,2,\ldots,K$. Denote the estimate by $\check{\gamma}^{(K)}(s,t)$;

     * 2.

Estimate $\sigma^{2}$ by
$\check{\sigma}_{K}^{2}(t)=\widetilde{r}^{(K)}(t)-\check{\gamma}^{(K)}(t,t)$,
where $\widetilde{r}^{(K)}(t)$ is obtained when estimating $\mu$ in step 2 of
(a), and note that it is assumed that the noise variance is a constant, we
estimate $\sigma^{2}$ by
$\check{\sigma}_{K}^{2}=\int\check{\sigma}_{K}^{2}(t)dt$;

     * 3.

Smooth
$\\{\widehat{C}_{Ki,j_{1},j_{2}}-\check{\gamma}^{(K)}(T_{Kij_{1}},T_{Kij_{2}})\\}^{2}$
by an online local linear smoother with the above $h_{\nu_{\gamma}}^{(K)}$ and
$\\{\eta_{\nu_{\gamma},j}^{(K)}\\}_{j=1}^{J}$ to obtain the estimate
$\check{{\mathrm{E}}}^{(K)}\Phi(s)^{2}\Phi(t)^{2}$.

Based on these,
$\widetilde{V}_{1}^{(K)}(s,t)=\check{{\mathrm{E}}}^{(K)}\Phi(s)^{2}\Phi(t)^{2}+\check{\sigma}_{K}^{2}\\{\check{\gamma}^{(K)}(s,s)+\check{\gamma}^{(K)}(t,t)\\}+\check{\sigma}_{K}^{4}-\check{\gamma}^{(K)}(s,t)^{2}$.

  3. (c)

Mean function estimation for dense data. Recall that for dense data,

$\displaystyle\Gamma_{\mu}=R(W)r(t)/f(t)+\gamma(t,t)$

Estimating $\gamma$ based on the raw covariance is infeasible as $\mu$ is not
available. As in (a), $r(t)$ is estimated by $\widetilde{r}^{(K)}(t)$, we can
make use of the dense measurements for each subject $(k,i)$, pre-smooth
$\\{Y_{kij}\\}_{j=1}^{m_{ki}}$ to obtain an estimate of
$X_{ki}(t)+\Phi_{ki}(t)$ which gives the estimate $\widetilde{\sigma}_{K}^{2}$
of noise variance $\sigma^{2}$ (Zhang and Chen, 2007) and estimate
$\gamma(t,t)$ by $\widetilde{r}^{(K)}(t)-\widetilde{\sigma}_{K}^{2}$:

     * 1.

Adopt the same method of (a) to compute $\widetilde{r}^{(K)}(t)$;

     * 2.

Pre-smooth to estimate $\sigma^{2}$. Specifically, for each $j$ at block $K$,
apply local linear smooth to $\\{Y_{Kij}:j=1,\ldots,m_{Ki}\\}$ with bandwidth
$h=O_{p}\big{(}\overline{m}_{\mu,K}^{-1/5}\big{)}$ to obtain smoothed
$\widehat{Y}_{Kij}$, then
$\check{{\varepsilon}}_{Kij}=Y_{Kij}-\widehat{Y}_{Kij}$. To de-bias, let
$\widehat{{\varepsilon}}_{Kij}=\check{{\varepsilon}}_{Kij}-\sum_{j=1}^{m_{Ki}}\check{{\varepsilon}}_{Kij}/m_{Ki}$
and estimate $\sigma^{2}$ by

$\widetilde{\sigma}_{K}^{2}=\frac{\sum_{i=1}^{n_{K}}\sum_{j=1}^{m_{Ki}}\widehat{{\varepsilon}}_{kij}^{2}}{S_{K,1}}+\frac{S_{K-1,1}}{S_{K,1}}\widetilde{\sigma}_{K-1}^{2},$

then estimate $\gamma(t,t)$ by
$\widetilde{r}^{(K)}(t)-\widetilde{\sigma}_{K}^{2}$.

  4. (d)

Covariance function estimation for dense data. Recall that

$\displaystyle\Gamma_{\gamma}=$
$\displaystyle\\{1+I(s=t)\\}\left\\{\frac{R(W)^{2}V_{1}(s,t)}{f(s)f(t)}+\frac{R(W)C_{1}}{C_{0}^{1/2}}\frac{f(s)V_{2}(t,s)+f(t)V_{2}(s,t)}{f(s)f(t)}\right\\}+C_{1}^{2}V_{3}(s,t),$

where $V_{1},V_{2},V_{3}$ are defined as in (4) of Section 4 and can be
written as
$V_{1}(s,t)={\mathrm{E}}\Phi(s)^{2}\Phi(t)^{2}+\sigma^{2}\\{\gamma(s,s)+\gamma(t,t)\\}+\sigma^{4}-\gamma(s,t)^{2}$,
$V_{2}(s,t)={\mathrm{E}}\Phi(s)^{2}\Phi(t)^{2}+\sigma^{2}\gamma(t,t)-\gamma(s,t)^{2}$,
$V_{3}(s,t)={\mathrm{E}}\Phi(s)^{2}\Phi(t)^{2}-\gamma(s,t)^{2}$. One can solve
$C_{1}$ by $C_{1}=\\{\nu_{\gamma}/(\alpha^{2}(W)\theta_{\gamma})\\}^{1/6}$
based on approximations of $V_{1},V_{2}$ and $V_{3}$. The estimation of
$V_{1}(s,t),V_{2}(s,t)$ and $V_{3}(s,t)$ consists of the three parts below:

     * 1.

Adopt the same method as in step 1 of (b) to compute
$\check{\gamma}^{(K)}(s,t)$;

     * 2.

Estimate $\sigma^{2}$ by $\widetilde{\sigma}_{K}^{2}$ in step 2 of (c);

     * 3.

Adopt the method in step 3 of (b) to compute $\check{\gamma}^{(K)}(s,t)$.

The above processes to estimate $\theta_{\mu},\theta_{\gamma}$ and
$\nu_{\mu},\nu_{\gamma}$ involve further bandwidth selection. Take the
standard local linear regression as an instance. According to Wand and Jones
(1994) and Fan and Yao (1998), the optimal bandwidths for estimating
$\theta_{\mu},\theta_{\gamma}$ and $\nu_{\mu},\nu_{\gamma}$ are

$h_{\theta_{\mu}}^{(K)}=G_{\mu}^{*}S_{K,1}^{-1/7},\
h_{\theta_{\gamma}}^{(K)}=G_{\gamma}^{*}S_{K,2}^{-1/8},\
h_{\nu_{\mu}}^{(K)}=R_{\mu}^{*}S_{K,1}^{-1/5},\
h_{\nu_{\gamma}}^{(K)}=R_{\gamma}^{*}S_{K,2}^{-1/6},$

where $G_{\mu}^{*},G_{\gamma}^{*},R_{\mu}^{*},R_{\gamma}^{*}$ involve unknown
quantities depending on $\mu$ an $\gamma$. As literature on classical
bandwidth selection has demonstrated, the convergences of
$\widetilde{h}_{\mu}^{(K)},\widetilde{h}_{\gamma}^{(K)}$ in (18) hold as long
as the bandwidths for estimating $\theta_{\mu},\theta_{\gamma},\nu_{\mu}$ and
$\nu_{\gamma}$ satisfy, for $G_{\mu},G_{\gamma},R_{\mu},R_{\gamma}=O_{p}(1)$,

$\displaystyle h_{\theta_{\mu}}^{(K)}=G_{\mu}S_{K,1}^{-1/7},\
h_{\theta_{\gamma}}^{(K)}=G_{\gamma}S_{K,2}^{-1/8},\
h_{\nu_{\mu}}^{(K)}=R_{\mu}S_{K,1}^{-1/5},\
h_{\nu_{\gamma}}^{(K)}=R_{\gamma}S_{K,2}^{-1/6}.$ (B.1)

Hence it is usually adequate to select appropriate constants to compute
$h_{\theta_{\mu}}^{(K)},h_{\theta_{\gamma}}^{(K)},h_{\nu_{\mu}}^{(K)},h_{\nu_{\gamma}}^{(K)}$.
Based on extensive numerical experiments, we recommend to set
$G_{\mu},R_{\mu}\in[0.5,1]$ and $G_{\gamma},R_{\gamma}\in[0.7,1]$. Using the
same argument of deriving (19), we use the following candidate bandwidths for
$h_{\theta}^{(K)}$ and $h_{\nu}^{(K)}$, for $j=1,2,\ldots,J$,

$\displaystyle\eta_{\theta_{\mu},j}^{(K)}=\\{(J-j+1)/J\\}^{1/7}h_{\theta_{\mu}}^{(K)},\
\eta_{\nu_{\mu},j}^{(K)}=\\{(J-j+1)/J\\}^{1/5}h_{\nu_{\mu}}^{(K)},$
$\displaystyle\eta_{\theta_{\gamma},j}^{(K)}=\\{(J-j+1)/J\\}^{1/8}h_{\theta_{\gamma}}^{(K)},\
\eta_{\nu_{\gamma},j}^{(K)}=\\{(J-j+1)/J\\}^{1/6}h_{\nu_{\gamma}}^{(K)}.$
(B.2)

So far we have delineated the whole process of online estimation for mean and
covariance functions in functional data.

## REFERENCES

  * (1)
  * Anava et al. (2013) Anava, O., Hazan, E., Mannor, S. and Shamir, O. (2013), ‘Online learning for time series prediction’, Journal of Machine Learning Research 30, 172–184.
  * Bacher et al. (2009) Bacher, P., Madsen, H. and Nielsen, H. A. (2009), ‘Online short-term solar power forecasting’, Solar Energy 83(10), 1772–1783.
  * Cai and Yuan (2011) Cai, T. T. and Yuan, M. (2011), ‘Optimal estimation of the mean function based on discretely sampled functional data: Phase transition’, The Annals of Statistics 39(5), 2330–2355.
  * Dekel et al. (2012) Dekel, O., Ran, G. B., Shamir, O. and Xiao, L. (2012), ‘Optimal distributed online prediction using mini-batches’, Journal of Machine Learning Research 13(1), 165–202.
  * Duchi and Singer (2009) Duchi, J. C. and Singer, Y. (2009), ‘Efficient online and batch learning using forward backward splitting’, Journal of Machine Learning Research 10(18), 2899–2934.
  * Fan and Yao (1998) Fan, J. and Yao, Q. (1998), ‘Efficient estimation of conditional variance functions in stochastic regression’, Biometrika 85(3), 645–660.
  * Hiraoka et al. (2000) Hiraoka, K., Yoshizawa, S., Hidai, K.-i., Hamahira, M., Mizoguchi, H. and Mishima, T. (2000), Convergence analysis of online linear discriminant analysis, in ‘International Joint Conference on Neural Networks’, Vol. 3, pp. 387–391.
  * James et al. (2000) James, G., Hastie, T. and Sugar, C. (2000), ‘Principal component models for sparse functional data’, Biomatrika 87(3), 587–602.
  * Kim et al. (2007) Kim, T.-K., Wong, S.-F., Stenger, B., Kittler, J. and Cipolla, R. (2007), Incremental linear discriminant analysis using sufficient spanning set approximations, in ‘2007 IEEE Conference on Computer Vision and Pattern Recognition’, pp. 1–8.
  * Kong et al. (2016) Kong, D., Xue, K., Yao, F. and Zhang, H. H. (2016), ‘Partially functional linear regression in high dimensions’, Biometrika 103(1), 147–159.
  * Kong and Xia (2019) Kong, E. and Xia, Y. (2019), ‘On the efficiency of online approach to nonparametric smoothing of big data’, Statistica Sinica 29(1), 185–201.
  * Kristan et al. (2010) Kristan, M., Skocaj, D. and Leonardis, A. (2010), ‘Online kernel density estimation for interactive learning’, Image and Vision Computing 28(7), 1106–1116.
  * Langford et al. (2009) Langford, J., Li, L. and Zhang, T. (2009), ‘Sparse online learning via truncated gradient’, Journal of Machine Learning Research 10(2), 777–801.
  * Li and Hsing (2010) Li, Y. and Hsing, T. (2010), ‘Uniform convergence rates for nonparametric regression and principal component analysis in functional/longitudinal data’, The Annals of Statistics 38(6), 3321–3351.
  * Lin and Xi (2011) Lin, N. and Xi, R. (2011), ‘Aggregated estimating equation estimation’, Statistics and Its Interface 4(1), 73–83.
  * Pang et al. (2005) Pang, S., Ozawa, S. and Kasabov, N. (2005), Incremental linear discriminant analysis for classification of data streams, in ‘IEEE transactions on Systems, Man and Cybernetics, part B (Cybernetics)’, Vol. 35, pp. 905–914.
  * Paul and Peng (2009) Paul, D. and Peng, J. (2009), ‘Consistency of restricted maximum likelihood estimators of principal components’, The Annals of Statistics 37(5), 1229–1271.
  * Rice and Silverman (1991) Rice, J. A. and Silverman, B. W. (1991), ‘Estimating the mean and covariance structure nonparametrically when the data are curves’, Journal of the royal statistical society series b-methodological 53(1), 233–243.
  * Richard et al. (2009) Richard, C., Bermudez, J. C. M. and Honeine, P. (2009), Online Prediction of Time Series Data With Kernels, Vol. 57, pp. 1058–1067.
  * Schifano et al. (2016) Schifano, E. D., Wu, J., Wang, C., Yan, J. and Chen, M. H. (2016), ‘Online updating of statistical inference in the big data setting’, Technometrics 58(3), 393–403.
  * Wand and Jones (1994) Wand, M. P. and Jones, M. C. (1994), Kernel Smoothing, Chapman and Hall.
  * Xiao (2010) Xiao, L. (2010), ‘Dual averaging method for regularized stochastic learning and online optimization’, Journal of Machine Learning Research 11(1), 2543–2596.
  * Yang et al. (2019) Yang, H., Pan, Z. and Tao, Q. (2019), ‘Online Learning for Time Series Prediction of AR Model with Missing Data’, Neural Processing Letters 50(3), 2247–2263.
  * Yao et al. (2005a) Yao, F., Muller, H. G. and Wang, J. L. (2005a), ‘Functional data analysis for sparse longitudinal data’, Journal of the American Statistical Association 100(June), 577–590.
  * Yao et al. (2005b) Yao, F., Muller, H. and Wang, J. (2005b), ‘Functional linear regression analysis for longitudinal data’, The Annals of Statistics 33(6), 2873–2903.
  * Zhang and Chen (2007) Zhang, J. and Chen, J. (2007), ‘Statistical inferences for functional data’, The Annals of Statistics 35(3), 1052–1079.
  * Zhang and Wang (2016) Zhang, X. and Wang, J.-L. (2016), ‘From sparse to dense functional data and beyond’, The Annals of Statistics 44(5), 2281–2321.
  * Zhou et al. (2003) Zhou, A., Cai, Z., Wei, L. and Qian, W. (2003), M-kernel merging: Towards density estimation over data streams, in ‘Eighth International Conference on Database Systems for Advanced Applications’, pp. 285–292.