1 引言

    局部多项式回归是非参数回归的一种方法，主要是由于Nadaraya−Watson估计方法（也叫Local-Constant）的加权是基于整个样本点，而且往往在边界上的估计效果并不理想。解决上述问题的办法就是用一个变动的函数取代局部固定的权重，即局部线性回归。局部线性回归就是在待估计点\(x\)的领域内用一个线性函数\(Y_i=\beta_0+\beta_1X_i,X_I\in[x-h,x+h]\)，来取代\(Y_i\)的平均，\(\beta_0,\beta_1\)是局部参数。

2 Nadaraya-Watson regression

    对于回归函数\(E\big(y_i|x_i'=x\big)=m(x)\), 使用非参数方法估计\(m(x)\).

\[ E[y_i|x'_i=x]=\int y f_{y|x}(y;x)dy=\int y\frac{f_{yx}(x,y)}{f_x(x)}dy \]

    使用核密度估计得到\(\hat{f}_{y x}(x,y),\hat{f}_{x}(x)\)

\[ \begin{array}{l}\hat{f}_x(x)=\frac{1}{nh^k}\sum_{i=1}^n K_1(\frac{x-\mathbf{x}_i}{h});\\ \hat{f}_{yx}(x,y)=\frac{1}{nh^{k+1}}\sum_{i=1}^n K_2\left(\frac{(y,x')-(\mathbf{y}_i,\mathbf{x}_i)}{h}\right).\end{array} \]

定义\(K_{h}(\cdot)=K(\cdot/h)/h\), 其中

\[ \frac{\int y\hat{f}(x,y;\mathbf{h})\mathrm{d}y}{\hat{f}_X(x;h_1)}=\frac{\int y\frac{1}{n}\sum_{i=1}^n K_{h_1}(x-X_i)K_{h_2}(y-Y_i)\mathrm{d}y}{\frac{1}{n}\sum_{i=1}^n K_{h_1}(x-X_i)} \\ =\frac{\frac{1}{n}\sum_{i=1}^nK_{h_1}(x-X_i)\int yK_{h_2}(y-Y_i)\mathrm{d}y}{\frac{1}{n}\sum_{i=1}^n K_{h_1}(x-X_i)} \\ =\frac{\frac{1}{n}\sum_{i=1}^n K_{h_1}(x-X_i)Y_i}{\frac{1}{n}\sum_{i=1}^n K_{h_1}(x-X_i)} \\ =\frac{\sum_{i=1}^n K_h(x-X_i)}{\sum_{i=1}^n K_h(x-X_i)}Y_i \\ =\sum\limits_{i=1}^n W_i(x)Y_i \]

    Nadaraya-Watson 估计量或者Local-Constant Estimator.为：

\[ \widehat{m}_h(x)=\frac{\sum_{i=1}^n K_h(X_i-x)Y_i}{\sum_{i=1}^n K_h(X_i-x)}. \]

    Nadaraya-Watson 估计量的风险为：

\[ \begin{aligned}R(\hat{f}_{n},f)=\frac{h^4}{4}\left\{\int x^2K(x)dx\right\}^2\int\left\{\ddot{f}(x)+2\dot{f}(x)\frac{\dot{g}(x)}{g(x)}\right\}^2dx\\ +\frac{\sigma^2\int K^2(x)dx}{nh}\int\frac{1}{g(x)}dx+o(nh^{-1})+o(h^4)\end{aligned} \quad \tag {1} \]

\(g\)是\(x\)的密度函数。

    让上式等于0并进行求导求解\(h\)得：最优带宽

\[ h_{\mathrm{opt}}=n^{-1 / 5}\left[\frac{\sigma^2 \int K^2(x) d x \int \frac{1}{g(x)} d x}{\left\{\int x^2 K(x) d x\right\}^2 \int\left\{\ddot{f}(x)+2 \dot{f}(x) \frac{\dot{g}(x)}{g(x)}\right\}^2 d x}\right]^{1 / 5} \]

Bias:

\[ h^2\left(\frac{1}{2} \ddot{f}(x)+\frac{\dot{f}(x) \dot{g}(x)}{g(x)}\right) \int u^2 K(u) d u+o\left(h^2\right) \]

方差为:

\[ \frac{\sigma^2(x)}{g(x) n h} \int K^2(u) d u+o\left(\frac{1}{n h}\right) \]

均方误差准则:

\[ \frac{1}{n h} \nu_0 \frac{\sigma^2(x)}{f_0(x)}+\frac{h^4}{4} g_0^{\prime \prime}(x)^2 \mu_2^2 \]

    Nadaraya-Watson 估计量等价于在每个点的局部通过求\(y\) 的均值来估计\(E(y|x)\)。使用”均值即常数回归”这一思路推广该估计方法：在\(x\)每个点的局部通过对进行多项式回归后取截距常数项作为\(E(y|x)\)的估计量，这样可以剔除局部的非线性因素。

3 local linear

3.1 Linear Regression with Weights

    对于Linear Regression with Weights来说，加权最小二乘法的目标是：

\[ r_w(\beta)=\sum_{i=1}^n w_i \varepsilon_i^2=\sum_{i=1}^n w_i\left(y_i-(X \beta)_i\right)^2 \]

写成矩阵形式:

\[ \begin{aligned} & r_w(\beta)=(y-X \beta)^{\top} W(y-X \beta) \\ & W=\left(\begin{array}{ccccc} w_1 & 0 & 0 & \cdots & 0 \\ 0 & w_2 & 0 & \cdots & 0 \\ 0 & 0 & w_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & w_n \end{array}\right) \end{aligned} \]

通过求导得：

\[ \hat{\beta}=\left(X^{\top} W X\right)^{-1} X^{\top} W y . \]

\[ \hat{\beta}_0+\hat{\beta}_1 \tilde{x}_1+\cdots+\hat{\beta}_p \tilde{x}_p=\tilde{x}^{\top} \hat{\beta} \]

where \(\tilde{x}=\left(1, \tilde{x}_1, \ldots, \tilde{x}_n\right)\).

3.2 Local Linear Regression

    为了描述 Local Linear Regression 的估计值 : 考虑 \(\mathrm{p}=1\) 是情形, 则 \(g(x)=g+x \beta\), 定义权函数 \(K_h(u)=h^{-r} K(u / h), w_i=K_h\left(x-x_i\right), \widehat{g}(x)\) 由下式给出:

\[ \min _{g, \beta} \sum_{i=1}^n\left(Y_i-g-\left(x-x_i\right)^{\prime} \beta\right)^2 K_h\left(x-x_i\right) . \]

\(\widehat{g}(x)\) 是加权最小二乘回归中的常数项, 回归形式为:

\[ \begin{gathered} Y=\left(\begin{array}{c} Y_1 \\ \vdots \\ Y_n \end{array}\right), X=\left(\begin{array}{cc} 1 & \left(x-x_1\right)^{\prime} \\ \vdots & \vdots \\ 1 & \left(x-x_n\right)^{\prime} \end{array}\right) \\ W=\operatorname{diag}\left(K_h\left(x-x_1\right), \ldots, K_h\left(x-x_n\right)\right) \end{gathered} \]

（1）参数估计为

\[ \begin{array}{l}\widehat{g}(x)=\widehat{\beta}_0\\=e_0^{\top}\hat{\beta}\\ =e_0^{\top}(X^{\top}WX)^{-1}X^{\top}Wy,\end{array} \]

\(e_0=(1,0,\ldots,0) \in \mathbb{R}^{p+1}\)，第一个位置为1，其他位置为0

（2）偏差、方差和MSE

    这个估计量是数据的局部线性拟合。对于距离\(x\)较远的观测值，它进行一个权重较小的回归。像核回归一样，这个估计量可以被解释\(Y_i\)为观测值的加权平均值，尽管权重有点复杂。

\[ \begin{gathered} S_0=\sum_{i=1}^n K_h\left(x-x_i\right), S_1=\sum_{i=1}^n K_h\left(x-x_i\right)\left(x-x_i\right), S_2=\sum_{i=1}^n K_h\left(x-x_i\right)\left(x-x_i\right)\left(x-x_i\right)^{\prime} \\ \hat{m}_0=\sum^n K_h\left(x-x_i\right) Y_i, \hat{m}_1=\sum^n K_h\left(x-x_i\right)\left(x-x_i\right) Y_i \end{gathered} \]

Bias

对于 bias, 考虑泰勒展开:

\[ g\left(x_i\right)=g_0(x)+g_0^{\prime}(x)\left(x_i-x\right)+\frac{1}{2} g_0^{\prime \prime}(x)\left(x_i-x\right)^2+\frac{1}{6} g_0^{\prime \prime \prime}\left(\bar{x}_i\right)\left(x_i-x\right)^3 . \]

Let \(r_i=g_0\left(x_i\right)-g_0(x)-\left[d g_0(x) / d x\right]\left(x_i-x\right)\).

用 \(\mathrm{X}\) 的形式表示:

\[ g=\left(g_0\left(x_1\right), \ldots, g_0\left(x_n\right)\right)^{\prime}=g_0(x) W e_1-g_0^{\prime}(x) W e_2+r \]

因为 \(e_1^{\prime} e_2=0\), 所以 bias 为:

\[ \begin{aligned} & e_1^{\prime}\left(X^{\prime} W X\right)^{-1} X^{\prime} W g-g_0(x)=e_1^{\prime}\left(X^{\prime} W X\right)^{-1} X^{\prime} W X e_1 g_0(x)-g_0(x) \\ & +e_1^{\prime}\left(X^{\prime} W X\right)^{-1} X^{\prime} W X e_2 g_0^{\prime}(x)+e_1^{\prime}\left(X^{\prime} W X\right)^{-1} X^{\prime} W r=e_1^{\prime}\left(X^{\prime} W X\right)^{-1} X^{\prime} W r \\ & e_1^{\prime}\left(X^{\prime} W X\right)^{-1} X^{\prime} W r=h^2 e_1^{\prime} H^{-1} \frac{\left(H^{-1} X^{\prime} W X H^{-1}\right)^{-1}}{n} \cdot \frac{h^{-2} H^{-1} X^{\prime} W r}{n} \\ & =\frac{h^2}{2} g_0^{\prime \prime}(x) e_1^{\prime}\left(\begin{array}{cc} 1 & 0 \\ 0 & \mu_2 \end{array}\right)^{-1}\left(\begin{array}{l} \mu_2 \\ \mu_3 \end{array}\right)=\frac{h^2}{2} g_0^{\prime \prime}(x) \mu_2 \text {. } \\ & \end{aligned} \]

同时, 这个估计量通过权重 \(K_h\left(x-x_i\right)\) 和回归量 \(x-x_i\) 依赖于 \(x\) 。所以计算 量会比较大。

方差:

\[ \begin{aligned} & e_1^{\prime}\left(X^{\prime} W X\right)^{-1} X^{\prime} W \Sigma W X\left(X^{\prime} W X\right)^{-1} e_1 \\ & =n^{-1} h^{-1} e_1^{\prime} H^{-1}\left(\frac{H^{-1} X^{\prime} W X H^{-1}}{n}\right)^{-1} \frac{h H^{-1} X^{\prime} W \Sigma W X H^{-1}}{n}\left(\frac{H^{-1} X^{\prime} W X H^{-1}}{n}\right)^{-1} H^{-1} e_1 \\ & =n^{-1} h^{-1}\left[\left(e_1^{\prime}\left[\begin{array}{cc} 1 & 0 \\ 0 & \mu_2 \end{array}\right]^{-1}\left[\begin{array}{cc} \nu_0 & \nu_1 \\ \nu_1 & \nu_2 \end{array}\right]\left[\begin{array}{cc} 1 & 0 \\ 0 & \mu_2 \end{array}\right]^{-1} e_1\right) \frac{\sigma^2(x)}{f(x)}+o_p(1)\right] \\ & \end{aligned} \]

MSE

\[ \frac{1}{n h} \nu_0 \frac{\sigma^2(x)}{f_0(x)}+\frac{h^4}{4} g_0^{\prime \prime}(x)^2 \mu_2^2 . \]

4 选择最优带宽

4.1 The leave-one-out cross-validation score:

\[ \mathrm{CV}(h)=\hat{R}(h)=\frac{1}{n} \sum_{i=1}^n\left\{Y_i-\hat{f}_{-i}\left(x_i\right)\right\}^2 \]

    \(\hat{f}_{-i}\left(x_i\right)\) 是去掉第 \(\mathrm{i}\) 个观测值得到的估计值。

    优点: 不受测试集合训练集划分方法的影响, 因为每一个数据都单独的做过 测试集。同时, 其用了 \(n-1\) 个数据训练模型, 也几乎用到了所有的数据, 保证 了模型的 bias 更小。不过 LOOCV 的缺点也很明显, 那就是计算量过于大,

4.2 The generalized cross-validation score:]

\[ \operatorname{GCV}(h)=\frac{1}{n} \sum_{i=1}^n\left\{\frac{Y_i-\hat{f}_n\left(x_i\right)}{1-\frac{1}{n} \operatorname{tr}(\mathbf{S})}\right\}^2 \]

4.3 AIC和BIC准则

R包-Loess

Loess局部加权多项式回归

    LOWESS最初由Cleveland 提出,后又被Cleveland&Devlin及其他许多人发展。在R中loess 函数是以lowess函数为基础的更复杂功能更强大的函数。主要思想为:在数据集合的每一点用低维多项式拟合数据点的一个子集,并估计该点附近自变量数据点所对应的因变量值,该多项式是用加权最小二乘法来拟合;离该点越远,权重越小,该点的回归函数值就是这个局部多项式来得到,而用于加权最小二乘回归的数据子集是由最近邻方法确定。

    最大优点：不需要事先设定一个函数来对所有数据拟合一个模型。并且可以对同一数据进行多次不同的拟合，先对某个变量进行拟合，再对另一变量进行拟合，以探索数据中可能存在的某种关系，这是普通的回归拟合无法做到的。

LOESS平滑方法

  1. 以\(x_0\)为中心确定一个区间，区间的宽度可以灵活掌握。具体来说，区间的宽度取决于\(q=fn\)。其中\(q\)是参与局部回归观察值的个数，\(f\)是参加局部回归观察值的个数占观察值个数的比例,\(n\)是观察值的个数。在实际应用中，往往先选定\(f\)值，再根据\(f\)和\(n\)确定\(q\)的取值，一般情况下\(f\)的取值在1/3到2/3之间。\(q\)与\(f\)的取值一般没有确定的准则。增大\(q\)值或\(f\)值，会导致平滑值平滑程度增加，对于数据中前在的细微变化模式则分辨率低，但噪声小，而对数据中大的变化模式的表现则比较好；小的\(q\)值或\(f\)值，曲线粗糙，分辨率高，但噪声大。没有一个标准的\(f\)值，比较明智的做法是不断的调试比较。

  2. 定义区间内所有点的权数，权数由权数函数来确定，比如立方加权函数\(weight = (1 - (dist/maxdist)^3)^3)\),dist为距离\(x\)的距离，maxdist为区间内距离\(x\)的最大距离。任一点\((x_0,y_0)\)的权数是权数函数曲线的高度。权数函数应包括以下三个方面特性：(1)加权函数上的点\((x_0,y_0)\)具有最大权数。(2)当\(x\)离开\(x_0\)时，权数逐渐减少。(3)加权函数以\(x_0\)为中心对称。

  3. 对区间内的散点拟合一条曲线\(y=f(x)\)。拟合的直线反映直线关系，接近\(x_0\)的点在直线的拟合中起到主要的作用，区间外的点它们的权数为零。

  4. \(x_0\)的平滑点就是\(x_0\)在拟合出来的直线上的拟合点\((y_0,f( x_0))\)。

  5. 对所有的点求出平滑点，将平滑点连接就得到Loess回归曲线。

loess(formula, data, weights, subset, na.action, model = FALSE,
       span = 0.75, enp.target, degree = 2,
       parametric = FALSE, drop.square = FALSE, normalize = TRUE,
       family = c("gaussian", "symmetric"),
       method = c("loess", "model.frame"),
       control = loess.control(...), ...)

formula是公式，比如y~x,可以输入1到4个变量;

data是放着变量的数据框，如果data为空，则在环境中寻找;

na.action指定对NA数据的处理，默认是getOption(“na.action”);

model是否返回模型框； span是alpha参数，可以控制平滑度,相当于上面所述的f，对于alpha小于1的时候，区间包含alpha的点，加权函数为立方加权，大于1时，使用所有的点最大距离为alpha^(1/p)，p 为解释变量;

anp.target，定义span的备选方法；

normalize，对多变量normalize到同一scale；

family,如果是gaussian则使用最小二乘法，如果是symmetric则使用双权函数进行再下降的M估计；

method，是适应模型或者仅仅提取模型框架；

control进一步更高级的控制，使用loess.control的参数；

    我们以cars数据集做为例子来看下使用方法。该数据中speed表示行驶速度，dist表示刹车距离。用loess来建立模型时重要的两个参数是span和degree，span表示数据子集的获取范围，取值越大则数据子集越多，曲线越为平滑。degree表示局部回归中的阶数，1表示线性回归，2表示二次回归，也可以取0，此时曲线退化为简单移动平均线。这里我们设span取0.4和0.8，从下图可见取值0.8的蓝色线条较为平滑。

plot(cars,pch=19)
model1=loess(dist~speed,data=cars,span=0.4)
lines(cars$speed,model1$fit,col='red',lty=2,lwd=2)
model2=loess(dist~speed,data=cars,span=0.8)
lines(cars$speed,model2$fit,col='blue',lty=2,lwd=2)
图1 Loess Function - Result
图1 Loess Function - Result

    当模型建立后，也可以类似线性回归那样进行预测和残差分析

x=5:25
predict(model2,data.frame(speed=x))
plot(model2$resid~model2$fit)
图二 Loess Function – Fit and Residual
图二 Loess Function – Fit and Residual

    R语言中另一个类似的函数是lowess，它在绘图上比较方便，但在功能上不如loess强大和灵活。

plot(cars,pch=19)
lines(lowess(cars),lty=2,lwd=2)
图三 lowess Function - Result
图三 lowess Function - Result

    LOESS的优势是并不需要确定具体的函数形式，而是让数据自己来说话，其缺点在于需要大量的数据和运算能力。LOESS作为一种平滑技术，其目的是为了探寻响应变量和预测变量之间的关系，所以LOESS更被看作一种数据探索方法，而不是作为最终的结论。

图四 lowess Function – Fit and Residual
图四 lowess Function – Fit and Residual

Reference

Avery, Matthew. Literature Review for Local Polynomial Regression. https://casual-inference.com/post/local-polynomial-smoothing.pdf. Accessed 7 June 2023.

Hughes, John. Local Polynomial Regression. 2022. https://johnhughes.org/bsta397/notes/Local_Polynomial_Regression.pdf. Accessed 7 June 2023.

“Locally Linear Regression.” MIT OpenCourseWare, http://ocw.mit.edu/courses/14-385-nonlinear-econometric-analysis-fall-2007/resources/local_lin_reg/. Accessed 7 June 2023.

The Local Polynomial Regression Estimator Deriving the Estimator. http://www.ma.man.ac.uk/~peterf/MATH38011/NPR%20local%20Linear%20Estimator.pdf. Accessed 7 June 2023.