关于局部加权回归（Locally Weighted Scatterplot Smoothing，LOWESS）、STL（Seasonal-Trend decomposition procedure b...

1. LOWESS（Locally Weighted Scatterplot Smoothing，局部加权回归）

0x1：lowess算法主要解决什么问题

1. 非线性回归拟合问题

LOWESS 通过取一定比例的局部数据，在这部分子集中拟合多项式回归曲线，这样我们便可以观察到数据在局部展现出来的局部规律和局部趋势（局部均值μ回归）。

同时，因为 LOWESS 的统计窗口缩小为局部窗口，因此拟合的回归曲线可以包含周期性，波动性的信息。

和传统的多元线性回归最大的不同是，通常的回归分析往往是根据全体数据建模，这样可以描述整体趋势，但现实生活中规律不总是（或者很少是）一条直线。

LOWESS将局部范围从左往右依次推进，最终一条连续的曲线就被计算出来了。显然，曲线的光滑程度与我们选取数据比例（局部窗口）有关：

• 局部统计窗口比例越少，拟合越不光滑，因为过于看重局部性质。但相对的，拟合的精确度更高
• 局部统计窗口比例越大，拟合越光滑，因为更看重整体性质。但拟合的精确度也相应降低

下面通过一个例子来形象说明"LOWESS"和"Linear regression"在局部信息捕获上的主要区别

数据文件是关于“物种数目与海拔高度”（中科院植物所-赖江山博士），数据集可以在这里下载。

用回归模型拟合的图示如下，本文采用R语言：

# 从counts.txt文件直接将数据读入R
x = read.csv("/Users/zhenghan/Downloads/counts.txt")
print(x)

par(las = 1, mar = c(4, 4, 0.1, 0.1), mgp = c(2.5,
                                              1, 0))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
abline(lm(counts ~ altitude, x), lwd = 2, col = "red")

如果仅仅从回归直线来看，似乎是海拔越高，则物种数目越多。如此推断下去，恐怕月球或火星上该物种最多。这显然是一个不是那么完美的线性回归推断。

改用lowess进行局部加权回归拟合，图示如下：

# 从counts.txt文件直接将数据读入R
x = read.csv("/Users/zhenghan/Downloads/counts.txt")
print(x)

par(las = 1, mar = c(4, 4, 0.1, 0.1))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
# 取不同的f参数值
for (i in seq(0.01, 1, length = 100)) {
  lines(lowess(x$altitude, x$counts, f = i), col = gray(i),
        lwd = 1.5)
  Sys.sleep(0.15)
}

上图中，曲线颜色越浅表示所取数据比例越大，即数据点越集中。不难看出白色的曲线几乎已呈直线状，而黑色的线则波动较大。

总体看来，图中大致有四处海拔上的物种数目偏离回归直线较严重：450 米、550 米、650 米和 700 米附近。这种偏离回归的现象恰好说明了实际问题的概率分布并不是严格遵循线性回归模型，而是其他的更复杂的非线性模型。

基于LOWESS的拟合结果，若研究者的问题是，多高海拔处的物种数最多？那么答案应该是在 650 米附近。 

为了确保我们用 LOWESS 方法得到的趋势是稳定的，我们可以进一步用 Bootstrap 的方法验证。因为 Bootstrap 方法是对原样本进行重抽样，根据抽得的不同样本可以得到不同的 LOWESS 曲线，最后我们把所有的曲线添加到图中，看所取样本不同是否会使得 LOWESS 有显著变化。

# 从counts.txt文件直接将数据读入R
x = read.csv("/Users/zhenghan/Downloads/counts.txt")
print(x)

set.seed(711) # 设定随机数种子，保证本图形可以重制
par(las = 1, mar = c(4, 4, 0.1, 0.1), mgp = c(2.5,
                                              1, 0))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
for (i in 1:400) {
  idx = sample(nrow(x), 300, TRUE) # 有放回抽取300个样本序号
  lines(lowess(x$altitude[idx], x$counts[idx]), col = rgb(0,
                                                          0, 0, 0.05), lwd = 1.5) # 用半透明颜色，避免线条重叠使得图形看不清
  Sys.sleep(0.05)
}

可以看出，虽然 700 米海拔附近物种数目减小的趋势并不明显了，这是因为这个海拔附近的观测样本量较少，在重抽样的时候不容易被抽到，因此在图中代表性不足，最后得到的拟合曲线分布稀疏。

但是经过 400 次重抽样并计算 LOWESS 曲线，刚才在第一幅图中观察到的趋势大致都还存在（因为默认取数据比例为 2/3，因此拟合曲线都比较光滑），所以整体回归拟合的结果是可信的。

2. 局部平滑问题

由于样本数据不可避免的随机波动问题（方差存在），我们获得的样本数据总是存在偏离均值的情况，为了更好地体现出数据中的统计规律性，我们需要对数据进行平滑处理。

举一个具体例子来说，在做数据平滑的时候，会有遇到有趋势或者季节性的数据，对于这样的数据，我们不能使用简单的均值正负3倍标准差以外做异常值剔除，需要考虑到趋势性等条件。使用局部加权回归，可以拟合一条趋势线，将该线作为基线，偏离基线距离较远的则是真正的异常值点。 

0x2：从KNN线性回归逐渐演进为LOWESS

上一章节讨论了线性回归在具体工程问题中面临的2个主要挑战，即：

• 周期性、波动性非线性数据回归问题
• 数据系统性噪点（非数据本身包含的方差，而是系统本身固有的方差）的平滑化问题

面对以上问题，学术界的陆续发表了很多研究paper，很好地解决了这些问题，笔者在这里尽量列举和梳理出整个理论体系的发展脉络。

1. 多元非线性回归模型

人们最早提出，用多元线性回归模型来拟合一元线性模型无法拟合的复杂样本数据。但是遇到了几个比较大的问题：

• 多元模型的模型函数需要人工设定，例如根据数据集绘出的图像，半观察半猜测要不要加个自变量的平方项，或者自变量的三角函数等等，作为新的回归算子加入，直到一些评价指标上拟合效果较好为止，这个过程就是非常繁重的特征工程工作
• 多元线性回归很容易陷入overfit或者underfit问题。

2. KNN线性回归

另一个解决问题的思路就是对数据进行局部最近邻处理，也即KNN线性回归。

其中：

• Nk(x)：为距离点x最近k个点xi组成的邻域集合（neighborhood set）
• yi：为KNN均值化后的回归结果
• ：为统计窗口内的回归结果，该统计窗口的所有点都用该点代替

但是，KNN线性回归有几个最大的缺点：

• 没有考虑到不同距离的邻近点应有不同的权重，不管是从统计上还是从领域经验角度来看，偏离均值越远的点重要性理应越低；
• 拟合的曲线不连续（discontinuous），拟合曲线的毛躁点还是很多，例如下图

 

关于KNN线性回归的问题的其他讨论，可以参阅这篇文章。

3. 引入kernel加权平滑的KNN线性回归

解决传统KNN线性回归问题的一个思考方向是：避免直接求均值这种“hard regression”的方式，改为寻找一种“soft regression”的方式。

我们在原始的KNN公式之上引入加权平滑：

权重因子函数K可以由很多种形式，一个大的原则就是：偏离均值越远的点，其贡献度相对越低。

比如，Epanechnikov二次kernel：

，

其中，λ为kernel的参数，称之为window width。对于KNN，只考虑最近的k个点影响；基于此，

，其中，x[k]为距离x0第k近的点，即窗口边界点。

经kernel加权平滑后，原始的KNN回归拟合的曲线为连续的了，如下图：

但是，这种kernel回归同样存在着边界（boundary）问题，如下图：

4. 基于最大似然估计的kernel加权平滑KNN回归

上一小节我们看到，传统kernel加权平滑的KNN回归，对于x序列的开始与结束区段的点，其左右邻域是不对称的，导致了平滑后的值偏大或偏小。

因此，需要对权值做再修正，修正的方法就是再次基于窗口内的数据点极大似然估计，使回归值x0进一步向均值线μ靠拢。

假定对统计窗口的回归量x0的估计值为：

定义目标函数：

令：

，

那么，目标函数可改写为：

求偏导，可得到：

利用矩阵运算，可直接得极大似然估计值为：

其中，

上述回归方法称之为LOWESS (Local Weighted regression)。

5. Robust LOWESS

我们沿着LOWESS的思路脉络继续讨论，截止目前为止，LOWESS似乎圆满地完成了所有任务，但是还有一个问题，就是样本数据中的异常点（outlier）问题。

我们知道，样本数据中的方差波动来自两个方面：

• 系统性误差：由于采样和测量手段引起的误差，不可避免
• 样本自身误差：由于样本数据自身存在的固有方差波动，导致的离群点

样本数据集的总体分布由以上两种方差波动共同作用产生，LOWESS通过kernel加权和极大似然回归的方式已经起到了较好的均值回归的作用，但是在极端的离群点（outlier）情况下，LOWESS的回归点x0依然可能偏离真值均值点μ过远。因此，人们又提出了Robust LOWESS方法。

Robust LOWESS是Cleveland在LOWESS基础上提出来的robust回归方法，能避免outlier对回归的影响。

在计算完估计值后，继续计算残差：

根据残差计算robustnest weight：

其中，s为残差绝对值序列的中位值（median），B函数为bisquare函数：

然后，用robustness weight乘以kernel weight作为的新weight。如此，便剔除了残差较大的异常点对于回归的影响，这本质上是数理统计中随机变量标准化的一种思想应用。

0x3：LOWESS对sin正弦周期数据的拟合效果测试

为了验证LOWESS的实际性能，我们来构造一组数据集验证其效果，该数据集的特性满足：

• 周期性：数据集不是单纯的线性回归关系，而是包含了周期波动性，我们基于sin函数生成y值
• 随机性：数据在不同的局部窗口中的密度分布是不同的，因此在不同局部窗口之间，均值并不是稳定连续变化的

lowess.py库函数

"""
lowess: Locally linear regression
==================================

Implementation of the LOWESS algorithm in n dimensions.

References
=========
[HTF] Hastie, Tibshirani and Friedman (2008). The Elements of Statistical
Learning; Chapter 6

[Cleveland79] Cleveland (1979). Robust Locally Weighted Regression and Smoothing
Scatterplots. J American Statistical Association, 74: 829-836.

"""
import numpy as np
import numpy.linalg as la


# Kernel functions:
def epanechnikov(xx, **kwargs):
    """
    The Epanechnikov kernel estimated for xx values at indices idx (zero
    elsewhere)

    Parameters
    ----------
    xx: float array
        Values of the function on which the kernel is computed. Typically,
        these are Euclidean distances from some point x0 (see do_kernel)

    Notes
    -----
    This is equation 6.4 in HTF chapter 6

    """
    l = kwargs.get('l', 1.0