3、LOESS回归基础：局部加权思想、权函数、平滑参数span

讲完为什么要做批次校正，我们得聊聊具体的工具了。

LOESS回归，全称是“局部加权散点图平滑”。说白了，它干的事就是——在你要预测的某个点附近，找一小片数据，做个简单的拟合。拟合完就扔掉，到了下一个点，再重新找一片邻居，重新拟合。

所以它根本没有全局公式，不像普通线性回归有个y = ax + b摆在那里。LOESS是“看着像什么就拟合什么”，我个人的理解是，它像一个特别灵活的手艺人，哪里不平补哪里。

3.1 局部加权思想

传统回归的毛病是什么？它假设所有数据点都是平等的。但LOESS不这么想。

举个我在项目中遇到的场景吧。测序数据里，往往有一批样本的某个基因表达量特别高，另一批特别低。你如果拿所有点一起拟合，结果就是高表达区被拉低，低表达区被拉高，两边都不讨好。LOESS的做法更聪明：在x₀这个位置做预测时，离x₀越近的点，权重越大；离得远，基本就不管了。这就是“局部加权”的实质。

• 邻居范围：由参数span控制，决定取多少比例的点参与拟合。
• 权重递减：不是一刀切，而是用一个函数平滑地降低远距离点的影响力。
• 拟合灵活：可以是直线，也可以是二次多项式，看数据复杂程度而定。

你想想看，这样做的好处是什么？它能很好地捕捉局部的数据趋势，不会因为远处几个离群点就把曲线带偏。嗯，这里要注意，局部加权的“局部”一词，其实暗含了一个假设：数据在局部区域是平滑变化的。如果你的数据像锯齿一样到处跳，那LOESS也救不了你。

3.2 权函数的选择

权重怎么给？这就得靠权函数了。我个人常用的有三种：tricube、gaussian，还有一个很少用的uniform。

先看tricube。这是LOESS的经典选项，也是推荐选项。它的公式长这样：

w(x) = (1 - |x|³)³   (当|x| ≤ 1时)
w(x) = 0              (当|x| > 1时)

说白了，就是用一个立方的立方来让权重从1平滑降到0。跟普通线性加权比，tricube降得更快、更圆滑。

# 实现tricube权函数，很简单
def tricube_weight(distance, bandwidth):
    scaled = distance / bandwidth
    if scaled > 1:
        return 0.0
    return (1 - scaled**3) ** 3

然后是gaussian权函数。它也用得挺多，尤其是当你希望权重永远不归零的时候。公式是：

w(x) = exp(- (x² / (2 * σ²)))

这里的σ通常被取成带宽的某个比例。跟tricube比，gaussian的尾巴拖得很长。即使离得很远的点，也有一点点权值。

我做过一个对比实验：相同的数据，用tricube和gaussian分别跑一遍LOESS，结果在边缘区域会有差异。gaussian会受远处样本的轻微影响，tricube则完全不管。哪个更好？看场景。如果你确定远处样本绝对不应该参与，就用tricube；如果你的数据噪音小，边缘点也可靠，gaussian也OK。

最后提一下uniform权函数，它的权重是常数，只要在带宽内，权值都一样。我个人觉得这有点浪费局部加权的优势，很少用。

3.3 平滑参数span

span是什么？它决定了一个点做拟合时会看周围多少个邻居。

具体来说，span是一个比例值，范围在0到1之间。span=0.3表示取该点附近30%的数据点来拟合；span=0.05就只取5%。

我习惯这样理解：span越小，对局部细节的拟合越“贴身”，但容易过拟合，把噪音也学进去了；span越大，曲线越平滑，但可能对真实结构欠拟合。

下面这张表能帮你快速理解span的影响：

项目里我通常会从0.3开始试。如果曲线太毛糙，就加大到0.5；如果发现拟合过头了（比如把噪音也当成趋势），就往回调到0.2。每次调整完之后，我都会拿残差图看一眼——残差如果呈现明显规律，说明span没选好。

为什么偏要用比例而不是固定点数？因为数据分布不均匀——高密度区可能需要较多点数，低密度区则少一些。用比例自动适配，省心。

3.4 动手试一试：简单实现

别光看理论，要动手。我用Python写过最简版的LOESS，核心代码也就几十行。下面是一个关键片段，展示如何用tricube权重和局部加权最小二乘做一个点：

import numpy as np

def loess_one_point(x_eval, x_data, y_data, span=0.3):
    # 计算所有点到当前点的距离
    distances = np.abs(x_data - x_eval)
    # 确定带宽：排序后取span比例的距离
    n = len(x_data)
    k = max(1, int(span * n))
    bandwidth = np.sort(distances)[k - 1]
    # 计算权重
    weights = np.where(distances <= bandwidth,
                       (1 - (distances / bandwidth)**3)**3,
                       0.0)
    # 加权最小二乘，这里用一阶多项式
    A = np.vstack([np.ones(n), x_data]).T
    W = np.diag(weights)
    # 解正规方程
    coeffs = np.linalg.lstsq(A.T @ W @ A, A.T @ W @ y_data, rcond=None)[0]
    return coeffs[0] + coeffs[1] * x_eval

这段代码里有一个关键点：每次预测一个新点，都要重新计算距离、带宽、权重、拟合。效率很低，但逻辑非常清晰。批量跑的时候可以向量化，我个人建议先理解单点逻辑，再考虑优化。