Kernel regression

Kernel regression is a class of nonparametric regression methods that estimates the conditional mean of a response variable by averaging observations according to their proximity to a target point. Unlike a fixed-dimensional linear model, it does not prescribe a global parametric form for the regression function. Its principal structural assumptions concern the smoothness of that function and the local distribution of the explanatory variables.

For observations ((X_i,Y_i)), (i=1,\ldots,n), the regression model is commonly written as

[ Y_i=m(X_i)+\varepsilon_i, \qquad \operatorname{E}(\varepsilon_i\mid X_i)=0, ]

where (m(x)=\operatorname{E}(Y\mid X=x)) is the unknown regression function. Kernel regression replaces (m(x)) by a locally weighted estimate whose behavior depends on a kernel function and a bandwidth.

Nadaraya–Watson estimator

The standard kernel regression estimator is the Nadaraya–Watson estimator,

[ \widehat m_h(x)

\frac{\displaystyle\sum_{i=1}^{n} K!\left(\frac{x-X_i}{h}\right)Y_i} {\displaystyle\sum_{i=1}^{n} K!\left(\frac{x-X_i}{h}\right)}, ]

where (K) is the kernel and (h>0) is the bandwidth. The denominator normalizes the local weights, so the estimate is a weighted average of the observed responses whenever the kernel is nonnegative. An observation receives greater weight when its explanatory value lies closer to the evaluation point under the distance scale induced by (h).

The estimator can also be written as a linear smoother,

[ \widehat m_h(x)=\sum_{i=1}^{n}W_i(x)Y_i, ]

with

[ W_i(x)= \frac{K((x-X_i)/h)} {\sum_{j=1}^{n}K((x-X_j)/h)}. ]

Although the estimator is linear in the observed responses, it is generally nonlinear as a function of the explanatory variables. Its normalization also distinguishes it from an unscaled convolution of the responses with a kernel.

The construction is closely related to kernel density estimation. If (\widehat f_h(x)) denotes a kernel estimate of the design density and (\widehat g_h(x)) denotes the corresponding kernel estimate of the response-weighted density, then the Nadaraya–Watson estimator has the ratio representation

[ \widehat m_h(x)=\frac{\widehat g_h(x)}{\widehat f_h(x)}. ]

This representation explains why the estimator becomes unstable in regions where the explanatory-variable density is small. The denominator contains limited local information in such regions, causing individual observations to exert comparatively large influence.

Historical development

The mathematical basis of kernel regression developed alongside the modern theory of density estimation. Murray Rosenblatt introduced a kernel density estimator in 1956, and Emanuel Parzen established a related formulation in 1962. In 1964, Elizbar Nadaraya and Geoffrey Watson independently described the normalized kernel estimator for conditional expectations that now bears their names.

Subsequent work connected kernel regression with local approximation and general smoothing theory. Charles Stone established broad consistency results for nonparametric regression procedures, while William Cleveland developed locally weighted regression into a systematic framework for statistical analysis. These developments placed kernel weighting within a wider class of methods based on fitting simple models over neighborhoods whose sizes contract as the sample grows.

Bias and variance

The bandwidth controls the estimator’s fundamental bias–variance balance. A small bandwidth concentrates weight in a narrow neighborhood, which reduces smoothing bias but leaves the estimate dependent on fewer observations. A large bandwidth incorporates more observations and lowers local sampling variability, while averaging over a wider portion of the regression function.

For a one-dimensional random design with density (f(x)), a symmetric second-order kernel, and a sufficiently smooth regression function, the leading interior bias of the Nadaraya–Watson estimator is

[ \operatorname{Bias}{\widehat m_h(x)} \approx \frac{h^2\mu_2(K)}{2} \left[ m''(x)+2m'(x)\frac{f'(x)}{f(x)} \right], ]

where

[ \mu_2(K)=\int u^2K(u),du. ]

The term involving (f'(x)) arises from the normalization by the estimated design density. Consequently, the local bias depends not only on the curvature of the regression function but also on changes in the concentration of explanatory-variable observations.

Under corresponding regularity conditions, the leading conditional variance is

[ \operatorname{Var}{\widehat m_h(x)\mid X_1,\ldots,X_n} \approx \frac{\sigma^2(x)R(K)} {nhf(x)}, ]

where

[ R(K)=\int K(u)^2,du ]

and (\sigma^2(x)=\operatorname{Var}(Y\mid X=x)). The variance therefore increases where the design density is low or the conditional response variance is high.

Balancing squared bias of order (h^4) against variance of order ((nh)^{-1}) gives an asymptotically optimal bandwidth of order (n^{-1/5}) for a twice-differentiable regression function in one dimension. The resulting pointwise mean squared error has order (n^{-4/5}). These rates change with the function’s smoothness, the dimension of the explanatory variable, and the form of the local estimator.

Local polynomial regression

Local polynomial regression generalizes the normalized local average by fitting a polynomial near each target point. At a location (x), a local polynomial of degree (p) minimizes

[ \sum_{i=1}^{n} K!\left(\frac{X_i-x}{h}\right) \left[ Y_i-\beta_0-\beta_1(X_i-x)-\cdots-\beta_p(X_i-x)^p \right]^2. ]

The fitted intercept (\widehat\beta_0) estimates (m(x)). When (p=0), the result is the Nadaraya–Watson estimator. When (p=1), the result is the local linear estimator.

Local linear regression removes the leading dependence of interior bias on the derivative of the design density. Under standard smoothness conditions, its principal interior bias is

[ \operatorname{Bias}{\widehat m_{\mathrm{LL}}(x)} \approx \frac{h^2\mu_2(K)}{2}m''(x). ]

It also responds differently near the boundary of the observed explanatory range. A local constant estimator has an asymmetric effective neighborhood at a boundary and consequently develops a first-order boundary bias. A local linear fit accounts for the direction and spacing of the available observations through its slope term, retaining second-order bias under ordinary boundary conditions.

Higher-degree local polynomials estimate derivatives as well as the regression function itself. Their equivalent weights need not remain nonnegative, because reproducing polynomial behavior requires moment constraints that differ from those of an ordinary weighted average.

Bandwidth selection and linear-smoother structure

Data-dependent bandwidth selection commonly uses estimates of prediction error. For leave-one-out cross-validation, the criterion has the form

[ \operatorname{CV}(h)

\frac{1}{n} \sum_{i=1}^{n} \left[ Y_i-\widehat m_{-i,h}(X_i) \right]^2, ]

where (\widehat m_{-i,h}) is calculated without the (i)-th observation. The selected bandwidth minimizes this criterion over a specified bandwidth family.

A 1983 analysis by You Watanabe expressed the leave-one-out calculation directly through the diagonal elements of the kernel smoother matrix. If the fitted values satisfy

[ \widehat{\mathbf m}=S_h\mathbf Y, ]

then the deleted residual obeys

[ Y_i-\widehat m_{-i,h}(X_i)

\frac{Y_i-\widehat m_h(X_i)} {1-(S_h)_{ii}}. ]

The identity reduces repeated refitting to a calculation based on one smoothing matrix for each bandwidth. It applies because the weights depend on the explanatory variables and bandwidth rather than on the response values.

The same linear-smoother representation underlies generalized cross-validation, which replaces individual diagonal corrections by an aggregate adjustment involving (\operatorname{tr}(S_h)). The trace acts as an effective measure of model complexity, although it is not generally an integer and does not correspond to a fixed number of regression coefficients.

Bandwidths may also be determined through asymptotic approximations to integrated mean squared error. Such criteria depend on unknown features of the design density, conditional variance, and regression derivatives. Plug-in methods estimate those quantities with preliminary smoothers before evaluating the bandwidth expression.

Multivariate form and dimensionality

For a (d)-dimensional explanatory variable, a bandwidth matrix (H) defines the local scale and orientation. The estimator becomes

[ \widehat m_H(x)

\frac{ \sum_{i=1}^{n} K!\left(H^{-1}(x-X_i)\right)Y_i }{ \sum_{i=1}^{n} K!\left(H^{-1}(x-X_i)\right) }. ]

A scalar bandwidth imposes the same smoothing scale in every direction, whereas a diagonal matrix permits coordinate-specific scales. A full positive-definite matrix additionally represents local orientation through correlations among coordinates.

The effective number of observations within a shrinking neighborhood decreases rapidly as dimension grows. Under an isotropic smoothness model with two derivatives, the conventional optimal bandwidth has order (n^{-1/(d+4)}), and the corresponding mean squared error has order (n^{-4/(d+4)}). This deterioration is an instance of the curse of dimensionality, reflecting the expanding data requirements of unrestricted multivariate smoothing.

Statistical interpretation

Kernel regression estimates a conditional expectation rather than a complete conditional distribution. It therefore describes the location of the response under squared-error loss but does not by itself determine conditional dispersion, asymmetry, or tail behavior. Related kernel methods estimate these features through nonparametric quantile regression, conditional density estimation, or separate smoothing of squared residuals.

The estimator is also sensitive to the geometry used to define proximity. Rescaling an explanatory variable changes kernel distances unless the bandwidth is transformed correspondingly. In multivariate applications, the bandwidth matrix therefore forms part of the statistical model’s effective geometry rather than serving only as a numerical tuning constant.

At points with sparse local support, normalized kernel estimates may exhibit high variance even when the kernel and bandwidth are fixed. Near gaps in the design distribution, observations on opposite sides of a low-density region may receive weight despite the absence of intermediate data. These properties follow from local averaging and are not represented by point estimates alone; inferential analysis commonly examines sampling distributions, confidence bands, or bootstrap approximations.