Universal kernel-type estimation of random fields ^††thanks: The study was supported by the program for fundamental scientific research of the Siberian Branch of the Russian Academy of Sciences, project FWNF-2022-0015.

AMS subject classification: 62G08.

1. Introduction

We study the following regression model:

$$Y_{i}=f(X_{i})+\xi_{i},\qquad i=1,\ldots,n,$$

(1)

where $f(t)$, $t:=(t^{(1)},\ldots,t^{(k)})\in\Theta\subset{\mathbb{R}}_{+}^{k}$, $k\geq 1$, is an unknown real-valued random function (a random field). We assume that $\Theta$ is a compact set, the random field $f(t)$ is continuous on $\Theta$ with probability 1, and the design $\{X_{i};\,i=1,\ldots,n\}$ consists of a collection of observed random vectors with unknown (generally speaking) distributions, not necessarily independent or identically distributed. The random design points $X_{i}$ may depend on $n$, i.e., a triangular array scheme for the design can be considered within this model. In particular, this scheme includes regression models with fixed design. Moreover, we do not require that the random field $f(t)$ be independent of the design $\{X_{i}\}$.

Next, we will assume that the unobservable random errors $\{\xi_{i}\}$ (a noise) form a martingale difference sequence, with

$$M_{p}:=\sup_{i}{\bf E}|\xi_{i}|^{p}<\infty\quad\mbox{for some}\,\,p>k\,\,\mbox{and}\,\,p\geq 2.$$

(2)

We also assume that $\{\xi_{i}\}$ are independent of the collection $\{X_{i}\}$ and the random field $f(\cdot)$. The noise $\{\xi_{i}\}$ may depend on $n$.

Our goal is to construct consistent in $C(\Theta)$ estimators for the random regression field $f(t)$ by the observations $\{(X_{i},Y_{i});\,i\leq n\}$ under minimal restrictions on the design points $\{X_{i}\}$, where $C(\Theta)$ denotes the space of all continuous functions on $\Theta$ with the uniform norm.

In the classical case of nonrandom $f(\cdot)$, the most popular estimating procedures are based on kernel-type estimators. We emphasize among them the Nadaraya–Watson, the Priestley-Chao, the Gasser-Müller estimators, the local polynomial estimators, and some of their modifications (see Härdle, 1990; Wand and Jones, 1995; Fan and Gijbels, 1996; Fan and Yao, 2003; Loader, 1999; Young, 2017; Müller, 1988). We do not aspire for providing a comprehensive review of this actively developing (especially in the last two decades) area of nonparametric estimation, and will focus only on publications representing certain methodological areas. We are primarily interested in conditions on the design elements. In this regard, a large number of publications in this area may be tentatively divided into the two groups: the studies dealing with fixed design $\{X_{i};\,i\leq n\}$ or with random one. In papers dealing with a random design, as a rule, the design consists of independent identically distributed random variables or stationary observations satisfying known forms of dependence, e.g., various types of mixing conditions, association, Markov or martingale properties, etc. Not attempting to present a comprehensive review, we may note the papers by Kulik and Lorek (2011), Kulik and Wichelhaus (2011), Roussas (1990, 1991), Györfi et. al. (2002), Masry (2005), Hansen (2008), Honda (2010), Laib and Louani (2010), Li et al. (2016), Hong and Linton (2016), Shen and Xie (2013), Jiang and Mack (2002), Linton and Jacho-Chavez (2010), Chu and Deng (2003) (see also the references in the papers). Besides, in the recent studies by Gao et al. (2015), Wang and Chan (2014), Chan and Wang (2014), Linton and Wang (2016), Wang and Phillips (2009a,b), Karlsen et al. (2007), the authors considered nonstationary design sequences under special forms of dependence (Markov chains, autoregressions, sums of moving averages, and so on).

In the case of fixed design, the vast majority of papers make certain assumptions on the regularity of design (see Zhou and Zhu, 2020; Benelmadani et al., 2020; Tang et al., 2018; Gu et al., 2007; Benhenni et al., 2010; Müller and Prewitt, 1993; Ahmad and Lin, 1984; Georgiev, 1988, 1990). In univariate models, the nonrandom design points $X_{i}$ are most often restricted by the formula $X_{i}=g(i/n)+o(1/n)$ with a function $g$ of bounded variation, where the error term $o(1/n)$ is uniform in $i=1,\ldots,n$. If $g$ is linear, then the design is equidistant. Another regularity condition in the univariate case is the relation $\max\nolimits_{i\leq n}(X_{i}-X_{i-1})=O(1/n)$, where the design elements are arranged in increasing order. In a number of recent studies, a more general condition $\max\nolimits_{i\leq n}(X_{i}-X_{i-1})\to 0$ can be found (e.g., see Yang and Yang, 2016; He, 2019; Wu et al., 2020). In several works, including those dealing with the so-called weighted estimators, certain conditions are imposed on the behavior of functions of design elements, but meaningful corresponding examples are limited to cases of regular design (e.g., see Zhang et al., 2019 ; Zhang et al., 2018; Liang and Jing, 2005; Roussas et al., 1992; Georgiev, 1988).

The problem of uniform approximation of the kernel-type estimators has been studied by many authors (e.g., see Einmahl and Mason, 2005; Hansen, 2008; Gu et al., 2007; Shen and Xie, 2013; Li et al., 2016; Liang and Jing, 2005; Wang and Chan, 2014; Chan and Wang, 2014; Gao et al., 2015 and the references therein).

In this paper, we study a class of kernel-type estimators, asymptotic properties of which do not depend on the design correlation structure. The design may be fixed (and not necessarily regularly spaced) or random (with not necessarily weakly dependent components). We present weighted least square estimators where the weights are chosen as the Lebesgue measures of the elements of a finite random partition of the regression function domain $\Theta$ such that every partition element corresponds to one design point. As a result, the proposed kernel estimators for the regression function are transformation of sums of weighted observations in a certain way with the structure of multiple Riemann integral sums, so that conceptually our approach is close to the methods of Priestley and Chao (1972) and of Mack and Müller (1988), who considered the cases of univariate fixed design and i.i.d. random design, respectively. Explicit upper bounds are obtained for the rate of uniform convergence of these estimators to the random regression field.

In contrast to the predecessors’ results, we do not impose any restrictions on the design correlation structure. We will consider the maximum cell diameter of the above-mentioned partition of $\Theta$ generated by the design elements, as the main characteristic of the design. Sufficient conditions for the consistency of the new estimators, as well as the windows’ widths will be derived in terms of that characteristic. The advantage of that characteristic over the classical weak dependence conditions is that the characteristic is insensitive to forms of correlation of the design elements. The main condition will be that the maximum cell diameter tends to zero in probability as the sample size grows. Note that such requirement is, in fact, necessary, since only when the design densely fills the regression function domain, it is possible to reconstruct the function more or less precisely.

Univariate versions of this estimation problem were studied in Borisov et al. (2021) and Linke et al. (2022) where the asymptotic analysis and simulations showed that the proposed estimators perform better than the Nadaraya–Watson ones in several cases. Note that the univariate case in Borisov et al. (2021) does not allow direct generalization to a multivariate case, since the weights were defined there as the spacings of the variational series generated by the design elements. Note also that the estimator in Borisov et al. (2021) is a particular univariate case of the estimators proposed in this paper, but not the only one. One of the univariate estimators studied here may be more accurate than the estimator in Borisov et al. (2021) (see Remark 3 below). Conditions on the design elements similar to those of this paper were used Linke and Borisov (2022), and in Linke (2023). The conditions provide uniform consistency of the estimators, but guarantee only pointwise consistency of the Nadarya–Watson ones. Besides, similar restrictions on the design elements were used before in Linke and Borisov (2017, 2018), and Linke (2019) in estimation of the parameters of several nonlinear regression models.

In this paper, we will assume that the unknown random regression function $f(t)$, $t\in\Theta$, is continuous with probability 1. Considering the general case of random regression function allows us to obtain results on estimating the mean function of a random regression process. In regard to estimating random regression functions, we may note the papers by Li and Hsing (2010), Hall et al. (2006), Zhou et al. (2018), Zhang and Wang (2016, 2018), Yao et al. (2005), Zhang and Chen (2007), Yao (2007), Lin and Wang (2022). In those papers, the mean and covariance functions of the random regression process $f(t)$ were estimated when, for independent noisy copies of the random process, each of the trajectories was observed in a certain subset of design elements (nonuniform random time grid). Estimation of mean and covariance of random processes is an actively developing area of nonparametric estimation, especially in the last couple of decades, is of independent interest, and plays an important role in subsequent analyses (e.g., see Hsing and Eubank, 2015; Li and Hsing, 2010; Zhang and Wang, 2016; Wang et al., 2016).

Estimation of random regression functions usually deals with either random or deterministic design. In the case of random design, it is usually assumed that the design elements are independent identically distributed (e.g., see Hall et al., 2006; Li and Hsing, 2010; Zhou et al., 2018; Yao , 2007; Yao et al., 2005; Zhang and Chen, 2007; Zhang and Wang, 2016, 2018; Lin and Wang, 2022). Some authors emphasized that their results can be extended to weakly dependent design (e.g., see Hall et al., 2006). For deterministic time grids, regularity conditions are often required (e.g., see Song et al., 2014; Hall et al., 2006). In regard to denseness of filling the regression function domain, the two types of design are distinguished in the literature: either the design is “sparse”, e.g., the number of design elements in each series is uniformly limited (Hall et al., 2006; Zhou et al., 2018; Li and Hsing, 2010), or the design is “dense” and the number of elements in a series increases with the sequential number of the series (Zhou et al., 2018; Li and Hsing, 2010). Uniform consistency of several estimators of the mean of random regression function was considered, for example, by Yao et al. (2005), Zhou et al. (2018), Li and Hsing (2010), Hsing and Eubank (2015), Zhang and Wang (2016).

In this paper, we consider one of the variants of estimation of the mean of a random regression function as an application of the main result. In the case of dense design, uniformly consistent estimators are constructed for the mean function, when the series-to-series-independent design is arbitrarily correlated inside each series. We require only that, in each series, the design elements form a refining partition of the domain of the random regression function. Our settings also include a general deterministic design situation, but we do not impose traditionally used regularity conditions. Thus, in the problem of estimating the mean function, as well as in the problem of estimating the function in model (1), we weaken traditional conditions on the design elements. Note that methodologies used for estimating the mean function for dense and for sparse data usually differ (e.g., see Wang et al., 2016). In the case of growing number of observations in each series, it is natural to preliminarily evaluate the trajectory of the random regression function in each series and then average the estimates over all series (e.g., see Hall et al., 2006). That is what we will do in this paper following this generally accepted approach. Universal estimates both for the mean and covariance functions of a random process in the case of sparse data, insensitive to the nature of the dependence of design elements, are proposed in Linke and Borisov (2024).

This paper has the following structure. Section 2 contains the main results on the rate of uniform convergence of the new estimators to the random regression function. In Section 3, we consider an application of the main results to the problem of estimating the mean and covariance function of a random regression field. In Section 4, the asymptotic normality of the new estimators is discussed. Section 5 contains several simulation examples. In Section 6, we discuss an example of assessing real data on earthquakes in Japan in 2012–2021. In Section 7, we summarize the main results of the paper. The proofs of the theorems and lemmas from Sections 2–4 are contained in Section 8.

2. Main assumptions and results

Without loss of generality we will assume that $d(\Theta\cup 0)\leq 1$, where

$$d(A):=\sup_{x,y\in A}\|x-y\|$$

is the diameter of a set $A$ and $\|\cdot\|$ is the supnorm in ${\mathbb{R}}^{k}$. In what follows, unless otherwise stated, all the limits will be taken as $n\to\infty$.

Our approach recalls a construction of multivariate Riemann integrals. To this end, we need the following condition on the design $\{X_{i}\}$.

$({\bf D})$ For each $n$, there exists a random partition of the set $\Theta$ into $n$ Borel-measurable subsets $\{\Delta_{i};\;i=1,\ldots,n\}$ such that $\delta_{n}:=\max_{i\leq n}d(\Delta_{i}\cup X_{i})\to 0$ in probability.

Condition (D) means that, for every $n$, the set $\{X_{i};\,i\leq n\}$ forms a $\delta_{n}$-net in the compact set $\Theta$. In particular, Condition (D) is satisfied if the design points $\{X_{i}\}$ are pairwise distinct, $X_{i}\in\Delta_{i}$ for all $i\leq n$, and $\max_{i\leq n}d(\Delta_{i})\to 0$ in probability.

In the case $\Theta=[0,1]^{k}$, a regularly spaced design satisfies Condition (D). Moreover, if $\{X_{i};\,i\geq 1\}$ is a stationary sequence satisfying an $\alpha$-mixing condition and $[0,1]^{k}$ is the support of the distribution of $X_{1}$, then Condition $(D)$ is fulfilled (see Remark 3 in Linke and Borisov, 2017). It is not hard to verify that, for i.i.d. design points with the probability density function of $X_{1}$ bounded away from zero on $[0,1]^{k}$, one can have $\delta_{n}=O\left(\frac{\log n}{n^{1/k}}\right)$ with probability 1. Notice that the dependence of random variables $\{X_{i}\}$ in Condition $(D)$ may be much stronger than that in these examples (see Linke and Borisov, 2017, 2018 and the example below).

Example. Let a sequence of bivariate random variables $\{{X}_{i};\,i\geq 1\}$ is defined by the relation

$${X}_{i}=\nu_{i}{U}_{1i}+(1-\nu_{i}){U}_{2i},$$

(3)

where the random vectors $\{{U}_{1i}\}$ and $\{{U}_{2i}\}$ are independent and uniformly distributed on the rectangles $[0,1/2]\times[0,1]$ and $[1/2,1]\times[0,1]$, respectively, while the sequence $\{\nu_{i}\}$ does not depend on $\{{U}_{1i}\}$ and $\{{U}_{2i}\}$ and consists of Bernoulli random variables with success probability $1/2$, i.e., the distribution of $X_{i}$ is the equi-weighted mixture of the two above-mentioned uniform distributions. The dependence between the random variables $\{\nu_{i}\}$ is defined by the equalities $\nu_{2i-1}=\nu_{1}$ and $\nu_{2i}=1-\nu_{1}$. In this case, the random variables $\{{X}_{i};\,i\geq 1\}$ in (3) form a stationary sequence of random variables uniformly distributed on the unit square $[0,1]\times[0,1]$, but, say, all known mixing conditions are not satisfied here because, for all natural $m$ and $n$,

$\displaystyle{\mathbb{P}}\big{(}{X}_{2m}\in[0,1/2]\times[0,1],\,{X}_{2n-1}\in[0,1/2]\times[0,1]\big{)}=0.$

On the other hand, it is easy to check that the stationary sequence $\{{X}_{i}\}$ satisfies the Glivenko–Cantelli theorem. This means that, for any fixed $h>0$,

$$\#\{i:\,\|{t}-{X}_{i}\|\leq h,\,1\leq i\leq n\}\sim 4h^{2}n$$

almost surely uniformly in $t$, where $\#$ denotes the standard counting measure. In other words, the sequence $\{{X}_{i}\}$ satisfies Condition $(D)$.

It is clear that, according to the scheme of this example, one can construct various sequences of dependent random variables uniformly distributed on $[0,1]\times[0,1]$, based on the choice of different sequences of the Bernoulli switches with the conditions $\nu_{j_{k}}=1$ and $\nu_{l_{k}}=0$ for infinitely many indices $\{j_{k}\}$ and $\{l_{k}\}$, respectively. In this case, Condition $(D)$ will also be satisfied. But the corresponding sequence $\{{X}_{i}\}$ (not necessarily stationary) may not satisfy the strong law of large numbers. For example, a similar situation occurs when $\nu_{j}=1-\nu_{1}$ for $j=2^{2k-1},\ldots,2^{2k}-1$ and $\nu_{j}=\nu_{1}$ for $j=2^{2k},\ldots,2^{2k+1}-1$, where $k=1,2,\ldots$ (i.e., we randomly choose one of the two rectangles $[0,1/2]\times[0,1]$ and $[1/2,1]\times[0,1]$, into which we randomly throw the first point, and then alternate the selection of one of the two rectangles by the following numbers of elements of the sequence: $1$, $2$, $2^{2}$, $2^{3}$, etc.). Indeed, we can introduce the notation $n_{k}=2^{2k}-1$, $\tilde{n}_{k}=2^{2k+1}-1$, $S_{m}=\sum\nolimits_{i=1}^{m}X_{i}^{(1)}$, with ${X}_{i}=(X_{i}^{(1)},X_{i}^{(2)})$, and note that, for all outcomes consisting the event $\{\nu_{1}=1\}$, one has

$$\frac{S_{n_{k}}}{n_{k}}=\frac{1}{n_{k}}\sum\limits_{i\in N_{1,k}}U_{1i}^{(1)}+\frac{1}{n_{k}}\sum\limits_{i\in N_{2,k}}U_{2i}^{(1)},$$

(4)

where ${U}_{ji}=(U_{ji}^{(1)},U_{ji}^{(2)})$, $j=1,2$; $N_{1,k}$ and $N_{2,k}$ are the collections of indices for which the observations $\{X_{i},i\leq n_{k}\}$ lie in the rectangles $[0,1/2]\times[0,1]$ or $[1/2,1]\times[0,1]$, respectively. It is easy to see that $\#(N_{1,k})=n_{k}/3$ and $\#(N_{2,k})=2\#(N_{1,k})$. Hence, ${S_{n_{k}}}/{n_{k}}\to{7}/{12}$ almost surely as $k\to\infty$ due to the strong law of large numbers for the sequences $\{U_{1i}^{(1)}\}$ and $\{u_{2i}^{(1)}\}$. On the other hand, for all elementary outcomes in the event $\{\nu_{1}=1\}$, as $k\to\infty$, we have with probability 1

$$\frac{S_{\tilde{n}_{k}}}{\tilde{n}_{k}}=\frac{1}{\tilde{n}_{k}}\sum\limits_{i\in\tilde{N}_{1,k}}U_{1i}^{(1)}+\frac{1}{\tilde{n}_{k}}\sum\limits_{i\in\tilde{N}_{2,k}}U_{2i}^{(1)}\to\frac{5}{12},$$

(5)

where $\tilde{N}_{1,k}$ and $\tilde{N}_{2,k}$ are the collections of indices for which the observations $\{{X}_{i},i\leq\tilde{n}_{k}\}$ lie the rectangles $[0,1/2]\times[0,1]$ or $[1/2,1]\times[0,1]$, respectively. In proving the convergence in (5) we took into account that $\#(\tilde{N}_{1,k})=(2^{2k+2}-1)/3$, $\#(\tilde{N}_{2,k})=2n_{k}/3$, i.e., $\#(\tilde{N}_{1,k})=2\#(\tilde{N}_{2,k})+1$.

Similar arguments are valid for elementary outcomes consisting the event $\{\nu_{1}=0\}$. $\hfill\Box$

In what follows, by $K(s)$, $s\in\mathbb{R}^{k}$, we will denote the kernel function. We assume that the kernel function is zero outside $[-1,1]^{k}$ and is a centrally symmetric probability density function, i.e., $K(s)\geq 0$, $K(s)=K(-s)$ for all $s\in[-1,1]^{k}$, and $\int_{[-1,1]^{k}}K(s)ds=1$. For example, we may consider product-kernels of the form

$$K(s)=\prod\limits_{j=1}^{k}K_{o}(s^{(j)}),$$

where $K_{o}(\cdot)$ is a univariate symmetric probability density function with support $[-1,1]$. We also assume that the function $K(s)$ satisfies the Lipschitz condition with constant $L\geq 1$:

$$|K(x)-K(y)|\leq L\left(|x^{(1)}-y^{(1)}|+\cdots+|x^{(k)}-y^{(k)}|\right)$$

for all $x=(x^{(1)},\dots,x^{(k)})$ and $y=(y^{(1)},\dots,y^{(k)})$, and put $K(y)=0$ for all $y$ such that $\|y\|>1$. Notice that, under these restrictions, $\sup_{s}K(s)\leq L$.

Put

$$K_{\varepsilon}(s):=\varepsilon^{-k}K(\varepsilon^{-1}s).$$

By $\theta_{\varepsilon}$ we denote a random vector with the density $K_{\varepsilon}(t)$, which is independent of the random variables $\{Y_{i}\}$.

Let $\Lambda(\cdot)$ denote the Lebesgue measure in ${\mathbb{R}}^{k}$. Introduce the following notation:

$\displaystyle f^{*}_{n,\varepsilon}(t):=\frac{\sum_{i=1}^{n}Y_{i}K_{\varepsilon}(t-X_{i})\Lambda(\Delta_{i})}{\sum_{i=1}^{n}K_{\varepsilon}(t-X_{i})\Lambda(\Delta_{i})},$

(6)

where $0/0=0$ by definition;

$\displaystyle J_{\varepsilon}(t):=\int\limits_{\Theta}K_{\varepsilon}(t-x)\,\Lambda(dx)\equiv{\bf P}(t-\theta_{\varepsilon}\in{\Theta}),\quad t\in\Theta;$

(7)

$$\omega_{f}(\varepsilon):=\sup_{x,y\in\Theta:\,\|x-y\|\leq\varepsilon}|f(x)-f(y)|.$$

Now, notice that

$${f^{*}_{n,\varepsilon}}(t)={\rm arg}\min\limits_{z\in\mathbb{R}}\sum\limits^{n}_{i=1}(Y_{i}-z)^{2}K_{\varepsilon}(t-X_{i})\Lambda(\Delta_{i}),$$

i.e., the estimators of the form (6) belong to the class of weighted least square estimators. Estimators (6) are also called local constant ones.

Finally, we will assume that there exist constants $\rho>0$ and $\varepsilon_{0}\in(0,1]$ such that

$$J_{\varepsilon}(t)\geq\rho\ \mbox{ for all }\ t\in\Theta\,\ \mbox{and positive}\,\,\varepsilon\leq\varepsilon_{0}.$$

(8)

So, some cases (for example, if $\Theta$ contains isolated points) are excluded from the scheme under consideration.

The main result of this paper is as follows.

Theorem 1. Let the conditions $(\rm D)$ and $(\ref{rho})$ hold. Then, for any fixed $\varepsilon\in(0,\varepsilon_{0}]$,

$\displaystyle\sup_{t\in\Theta}|f^{*}_{n,\varepsilon}(t)-f(t)|\leq\omega_{f}(\varepsilon)+\zeta_{n}(\varepsilon)$

(9)

with probability $1$, where $\zeta_{n}(\varepsilon)$ is a positive random variable such that

	$\displaystyle{\bf P}\left(\zeta_{n}(\varepsilon)>y\right)$	$\displaystyle\leq$	$\displaystyle G(k,p)\,\rho^{-p}\,M_{p}\,L^{p/2}\,y^{-p}\,\varepsilon^{-k(p/2+1)}\,{\bf E}(\delta_{n}^{kp/2})\,$		(10)
			$\displaystyle+\ {\bf P}(\delta_{n}>\varepsilon\min\{1,\,\rho(k2^{k+1}L)^{-1}\}),$

where

$$0<G(k,p)<(p-1)^{p/2}2^{p(k+(3/2))}\,\left(1+\frac{k}{2^{(p-k)/(p+1)}-1}\right)^{p+1}.$$

In what follows, we will denote by $O_{p}(\eta_{n})$ some univariate random variables $\zeta_{n}$ such that, for all $y>0$,

$$\limsup_{n\to\infty}{\bf P}(|\zeta_{n}|/\eta_{n}>y)\leq\beta(y),$$

where $\{\eta_{n}\}$ is a sequence of positive nonrandom numbers, $\lim_{y\to\infty}\beta(y)=0$, and the function $\beta(y)$ does not depend on $n$.

Taking Remark 1 into account one can obtain the following two assertions as consequences of Theorem $1$.

Corollary 1. Let $\cal C$ be a set of nonrandom equicontinuous functions from the function space $C[0,1]^{k}$ $($for example, a precompact subset of $C[0,1]^{k}$$)$. Then, under Condition $(D)$,

$$\gamma_{n}({\cal C}):=\sup_{f\in\cal C}\,\sup_{t\in[0,1]^{k}}|f^{*}_{n,\tilde{\varepsilon}(n)}(t)-f(t)|\stackrel{{\scriptstyle p}}{{\to}}0,$$

where $\tilde{\varepsilon}(n)$ is a solution to equation $(\ref{opt})$ in which the modulus of continuity $\omega_{f}(\varepsilon)$ is replaced with the universal modulus $\omega_{\cal C}(\varepsilon):=\sup_{f\in\cal C}\omega_{f}(\varepsilon)$. In this case, $\gamma_{n}({\cal C})=O_{p}(\omega_{\cal C}(\tilde{\varepsilon}(n)))$.

Corollary 2. If the modulus of continuity of the regression random field $f(t)$, $t\in[0,1]^{k}$, in Model $(\ref{f2})$ meets the condition $\omega_{f}(\varepsilon)\leq\zeta\tau(\varepsilon)$ a.s., where $\zeta>0$ is a proper random variable and $\tau(\varepsilon)$ is a positive continuous nonrandom function, with $\tau(\varepsilon)\to 0$ as $\varepsilon\to 0$, then, under Condition $(D)$,

$\displaystyle\sup_{t\in[0,1]^{k}}|f^{*}_{n,\hat{\varepsilon}(n)}(t)-f(t)|\stackrel{{\scriptstyle p}}{{\to}}0,$

(12)

where $\hat{\varepsilon}(n)$ is a solution to equation $(\ref{opt})$ in which the modulus of continuity $\omega_{f}(\varepsilon)$ is replaced with $\tau(\varepsilon)$.

Example 2. Let $\Theta=[0,1]^{k}$, $\delta_{n}\leq\nu n^{-1/k}$, with ${\bf E}\nu^{kp/2}<\infty$, and $\omega_{f}(\varepsilon)\leq\zeta\varepsilon^{\alpha}$, $\alpha\in(0,1]$, where $\zeta$ is a proper random variable. Then $\varepsilon(n)=O\left(n^{-\frac{1}{2k(1/p+1/2)+\alpha}}\right)$ and

$$\sup_{t\in[0,1]^{k}}|f^{*}_{n,\varepsilon}(t)-f(t)|=O_{p}\left(n^{-\frac{\alpha}{2k(1/p+1/2)+\alpha}}\right).$$

In particular, in the one-dimensional case, if $f(\cdot)=W(\cdot)$ is a Wiener process on $[0,1]$, and the i.i.d. random variables $\xi_{i}$ are centered Gaussian, then by Lévy’s modulus of continuity theorem, for any arbitrarily small $\nu>0$, we have

$$\sup_{t\in[0,1]}|f^{*}_{n,\varepsilon}(t)-f(t)|=O_{p}(n^{-1/3+\nu}).$$

Here we put $k=1$, $\alpha=1/2-\nu_{1}$, and $1/p<\nu_{2}$, with arbitrarily small positive $\nu_{1}$ and $\nu_{2}$.

3. Application to estimating the mean and covariance functions
of a random regression function

In this section, as an application of Theorem 1, we will construct a consistent estimator for the mean function of the random regression function in Model $(\ref{f2})$. We consider the following multivariate statement of the problem of estimating the mean function of an a.s. continuous random regression stochastic process. Consider $N$ independent copies of Model $(\ref{f2})$:

$$Y_{i,j}=f_{j}(X_{i,j})+\xi_{i,j},\qquad i=1,\ldots,n,\,\,\,\,j=1,\ldots,N,$$

(16)

where $f(t),f_{1}(t),\ldots,f_{N}(t)$, $t\in[0,1]^{k}$, are i.i.d. unknown a.s. continuous stochastic processes, and, for every $j$, the collection $\{\xi_{i,j};\,i\leq n\}$ satisfies condition (2). Here and in what follows, the subscript $j$ denotes the sequential number of such a copy. Introduce the notation

$$\overline{f^{*}}_{N,n,\varepsilon}(t):=\frac{1}{N}\sum\limits_{j=1}^{N}f^{*}_{n,\varepsilon,j}(t).$$

Theorem 2. Let Condition $(D)$ for Model $(\ref{f2})$ be fulfilled and

$${\bf E}\sup_{t\in[0,1]^{k}}|f(t)|<\infty.$$

(17)

Besides, let a sequences $\varepsilon\equiv\varepsilon_{n}\to 0$ and a sequence of naturals $N\equiv N_{n}\to\infty$ satisfy the conditions

$$\varepsilon^{-k(p/2+1)}\,{\bf E}(\delta_{n}^{kp/2})\to 0\,\,\,\mbox{and}\,\,\,N{\bf P}(\delta_{n}>\varepsilon\min\{1,\,\rho(k2^{k+1}L)^{-1}\})\to 0.$$

(18)

Then

$$\sup\limits_{t\in[0,1]^{k}}|\overline{f^{*}}_{N,n,\varepsilon}(t)-{\bf E}f(t)|\stackrel{{\scriptstyle p}}{{\to}}0.$$

(19)

In this section, we discuss sufficient conditions for asymptotic normality of the estimators $f^{*}_{n,\varepsilon}(t)$. Denote by ${\cal F}_{0}$ the trivial $\sigma$-field, and by ${\cal F}_{j}$ the $\sigma$-field generated by the collection $\{\xi_{1},\ldots,\xi_{j}\}$, by the design points, and by the regression random field.

Theorem 3. Let the design $\{X_{i}\}$ do not depend on $n$. Under Condition $(D)$, assume that, for some $t\in\Theta$ and a sequence $\varepsilon\equiv\varepsilon_{n}$,

$$h_{n}:=\frac{\max_{j\leq n}\big{(}K_{\varepsilon}\left(t-X_{j}\right)\Lambda(\Delta_{j})\big{)}^{2}}{\sum_{j=1}^{n}\big{(}K_{\varepsilon}\left(t-X_{j}\right)\Lambda(\Delta_{j})\big{)}^{2}}\stackrel{{\scriptstyle p}}{{\to}}0,$$

(20)

$${\bf E}(\xi_{j}^{2}\ |\ {\cal F}_{j-1})=\sigma^{2}\ \mbox{ a.s. for all }j,$$

$$\max_{j}{\bf E}\left(\xi_{j}^{2}\,{\bf 1}(\xi_{j}^{2}>a/h_{n})\ |\ {\cal F}_{j-1}\right)\stackrel{{\scriptstyle p}}{{\to}}0\ \mbox{ for all }a>0.$$

Then

$$B^{-1}_{n,\varepsilon}(t)\left(f^{*}_{n,\varepsilon}(t)-f(t)-r_{n,\varepsilon}(t)\right)\stackrel{{\scriptstyle d}}{{\to}}N(0,\sigma^{2}),$$

where $N(0,\sigma^{2})$ is a centered Gaussian random variable with variance $\sigma^{2}$,

$$B^{2}_{n,\varepsilon}(t):=J^{-2}_{n,\varepsilon}(t)\sum\limits_{i=1}^{n}\big{(}K_{\varepsilon}(t-X_{i})\Lambda(\Delta_{i})\big{)}^{2},$$

$$r_{n,\varepsilon}(t):=J^{-1}_{n,\varepsilon}(t)\sum\limits_{i=1}^{n}(f(X_{i})-f(t))K_{\varepsilon}\left(t-X_{i}\right)\Lambda(\Delta_{i}),$$

$$J_{n,\varepsilon}(t):=\sum\limits_{i=1}^{n}K_{\varepsilon}(t-X_{i})\Lambda(\Delta_{i}).$$

The theorem is a direct consequence of Corollary 3.1 in Hall and Heyde (1980).

In this section, we present simulations comparing the estimator $f^{*}_{n,\varepsilon}(t)$ defined in (6) with the Nadaraya–Watson estimator

$$\hat{f}_{n,\varepsilon}(t):=\frac{\sum_{i=1}^{n}Y_{i}K_{\varepsilon}(t-X_{i})}{\sum_{i=1}^{n}K_{\varepsilon}(t-X_{i})}$$

in the 2-dimensional case. For this estimator, we will assume $0/0=0$, like that was done for the estimator $(\ref{main})$.

The elements of the design space $\Theta=[-1,1]\times[-1,1]$ will be denoted by $(x,y)$. The following two algorithms were used to partition the space $\Theta$ into the sets $\Delta_{i}$.

The first algorithm is the Voronoi partitioning. For each $i$, the set $\Delta_{i}$ is the Voronoi cell corresponding to $X_{i}$, i.e. the set of all points of $\Theta$ that lie closer to $X_{i}$ than to any other design point. The deldir R package was employed for calculation of the squares of the cells.

The second algorithm is recursive partitioning by coordinate-wise medians. First, we divide $\Theta$ into the two rectangles by the line $t^{(1)}={\rm median}\{X_{1}^{(1)},\dots,X_{n}^{(1)}\},$ where the median is the midpoint of the interval $\left(X_{\lfloor n/2\rfloor}^{(1)},\ X_{\lfloor n/2\rfloor+1}^{(1)}\right)$ when all the points are sorted in increasing order with respect to the first coordinate. Then each of the two rectangles is divided recursively. If, at some step, a rectangle contains two or more design points then it is divided into the two parts: If the rectangle’s width is greater than its height, then the rectangle is divided by the line $t^{(1)}={\rm median}\{X_{l}^{(1)}:l\in B\},$ where $B$ is the set of indices of the design points falling into the rectangle; otherwise it is divided by the line $t^{(2)}={\rm median}\{X_{l}^{(2)}:l\in B\}$. As soon as there is only one design point $X_{i}$ in a rectangle, the rectangle is put to be $\Delta_{i}$.

Results of partitioning $\Theta$ into cells for a collection of 100 points by the both algorithms are displayed in Fig. 1.

In the simulation examples below, we used the tricubic kernel

$$K(x,y)=\frac{440}{162\pi}\max\left\{0,\left(1-\sqrt{x^{2}+y^{2}}^{\,3}\right)^{3}\right\}.$$

In each example, 1000 simulation runs were performed. In each of the simulation runs, 5000 design points were generated and randomly divided into the training (80%) and validation (20%) sets. For the design points $X_{i}$, the observations $Y_{i}$ were generated with i.i.d. Gaussian noise with standard deviation $\sigma=0.5$. For each of the tested algorithms, on the training set, the optimal $\varepsilon$ was calculated by 10-fold cross-validation minimizing the average of mean-square errors. The $\varepsilon$ was selected from 20 values located on the logarithmic grid from 0.01 to 0.5. The random partitioning for the cross-validation was the same for all the tested algorithms.

Then, for each of the algorithms, the model, trained on the training set with the chosen $\varepsilon$, was used to compute the mean-square error (MSE) for the observations of the validation set:

$${\rm MSE}=\frac{1}{m}\sum_{j}(f^{*}_{\varepsilon}(X_{j})-Y_{j})^{2},$$

where the sum is taken over the validation set, and $m$ is the size of the set. Besides, that model was employed to compute the maximal absolute error (MaxE) for the true values of the target function $f$ on the $100\times 100$ uniform lattice on $\Theta$:

$${\rm MaxE}=\max_{j}|f^{*}_{\varepsilon}(\gamma_{j})-f(\gamma_{j})|,$$

where the maximum is computed for the elements $\gamma_{j}$ of the $100\times 100$ lattice covering $\Theta$.

The algorithms that were compared will be denoted by NW (Nadaraya-Watson), ULCV (Universal Local Constant estimator (6) with Voronoi partitioning), and ULCM (Universal Local Constant estimator (6) with coordinate-wise Medians partitioning).

The results of the simulation runs are presented as median (1-st quartile, 3-rd quartile) and are compared between the estimators with the paired Wilcoxon test.

In the examples below, we intentionally chose the densities of the design points with high nonuniformity in order to demonstrate possible advantages of the new estimator.

In this example, we approximate the nonrandom regression function

$$f(x,y)=\frac{5}{1+e^{-20x}}-2y^{3}.$$

The design points were generated in a way similar to that in Example in Sec. 1. First, we choose the left rectangle $[-1,0)\times[-1,1]$ or the right rectangle $[0,1]\times[-1,1]$ with equal probabilities and draw $X_{1}$ uniformly distributed in the chosen rectangle. Then we draw $10$ design points uniformly distributed in the other rectangle. Then we draw $10^{2}$ design points uniformly distributed in the rectangle where $X_{1}$ lies. Then we draw $10^{3}$ design points uniformly distributed in the rectangle where $X_{2}$ lies, and so on. In other words, we alternate the rectangle after $1$, $10$, $10^{2}$, $10^{3}$, … draws. One draw of the $5000$ design points is depicted in Fig. 2. The estimated function $f(x,y)$ and a computed ULCV estimate are depicted in Fig. 3.

The results are presented in Fig. 2. The ULCV estimator appeared to perform best among the three considered ones both for MSE and MaxE accuracy measures. In particular, the ULCV estimator was better than the NW one: MSE 0.2661 (0.2584, 0.2742) vs. 0.2734 (0.2650, 0.2819), $p<0.0001$; MaxE 0.7878 (0.7013, 0.9230) vs. 1.1998 (1.0911, 1.3250), $p<0.0001$. In this example, the ULCM estimator was better than the NW one as well.

In this example, we approximate the nonrandom regression function

$$f(x,y)=\sin\left(10\sqrt{x^{2}+y^{2}}\right)/\sqrt{x^{2}+y^{2}}.$$

(21)

The i.i.d. design points were generated with independent polar coordinates $(\rho,\varphi)$, where $\rho$ was drawn with the density proportional to $r^{2}(2-r)^{1/10}$, $0\leq r\leq 2$, and $\varphi$ was uniformly distributed on $[0,2\pi)$. The distribution of the design points was restricted on $\Theta=[-1,1]\times[-1,1]$, i.e., the design points that did not fall into $\Theta$ were excluded, keeping the total number of collected points equal to 5000, as in the other simulation examples. One draw of the design points is depicted in Fig. 4. The estimated function $f(x,y)$ and a computed ULCV estimate are depicted in Fig. 5.