ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

线性回归实战:从最小二乘法到sklearn模型评估与避坑指南

线性回归实战:从最小二乘法到sklearn模型评估与避坑指南

先说一个我辅导新人时常遇到的问题。很多人拿到数据后的第一反应就是打开 sklearn,调用 LinearRegression,fit 完看一眼 R²,数字不错就觉得任务结束了。但换到真实业务数据,结果往往一言难尽:不是系数符号跟业务经验相反,就是换几组数据结果抖得厉害,又或者测试集表现和训练集完全对不上。这时候再回头看公式、补原理,成本已经很高了。

所以这篇内容我想换一种讲法:先用一个小例子把线性回归在解决的问题讲清楚,再跟着数据用手写最小二乘法把系数算出来,最后才用 sklearn 封装好的工具落地。整个过程会尽量保持“说人话”的风格,也把容易踩的坑提前摆出来。不管你是刚接触机器学习的初学者,还是已经写过不少模型但想回头补基础的人,这篇都可以当成一份完整的实操笔记来用。

1. 先想清楚:线性回归到底在拟合什么

1.1 用一个小例子找到模型直觉

假设你手里有一家小店的历史记录,想根据每天投放的广告曝光量和上架商品数量来预测当天的营业额。你手上有一批过去几十天的数据,每一条记录都包含三列:广告曝光量、上架商品数、当天实际营业额。

线性回归的思路很简单:给两个输入变量分别配上一个权重,再加上一个基础值,然后得到一个预测营业额的公式。比如预测值等于基础值 + 权重1 × 广告曝光量 + 权重2 × 上架商品数。训练模型的过程,就是在找一组“最合适”的权重和基础值,让公式算出来的预测营业额尽可能贴近真实营业额。

把这套逻辑推广到多个变量的场景,它就变成了数学里的一条超平面。线性回归这个名词里的“线性”指的不是说数据本身一定是直线关系,而是说我们假设输出和输入之间的关系可以被这些权重的线性组合描述。你可以把它理解成在平面上画一条最“贴近”所有点的直线,只不过当输入变多时,这条线变成了高维空间里的一个平面。

1.2 从数学符号到损失函数

用记号来表达,假设我们有 n 条样本、p 个特征,那么每条样本可以写成:

y_i = β0 + β1 * x_i1 + β2 * x_i2 + ... + βp * x_ip + ε_i

这里的 β0 是截距,β1 到 βp 是每个特征对应的系数,ε_i 是模型没有解释掉的误差。写成矩阵形式会更紧凑:

y = Xβ + ε

X 是一个 n 行 p+1 列的矩阵,第一列全是 1,用来对应截距项。如果只考虑线性回归本身,我们不需要对 ε 的分布做太严格的假设,但通常默认它的均值是 0。

那么问题来了:怎么判断一组 β 好不好?最自然的办法是计算每条样本的预测值和真实值之间的差距,然后把所有差距加起来。为了防止正负抵消,常见的做法是使用平方误差,于是就有了目标函数:

L(β) = Σ(y_i - X_iβ)^2

所有线性回归的核心,本质上就是找到一组让这个平方误差和最小的 β,也就是最小二乘估计。为什么用平方而不是直接用绝对值?一方面是因为平方函数处处可导,求导方便;另一方面是如果我们额外假设误差服从正态分布,最小化平方误差等价于最大化最大似然估计。虽然实际数据很难完美满足正态假设,但这个选择在大多数场景下都足够稳定。

1.3 简单但容易被忽略的假设

很多人学完公式直接就开始跑数据,却忽略了线性回归有一组隐含假设。不是说数据不满足假设就一定不能跑,而是你得知道自己正在接受哪些风险。

第一个假设是输入和输出之间的关系大致可以用线性来描述。如果真实关系是曲线,线性回归也会硬拟合出一条直线,但残差会非常难看。第二个是误差的均值要为零,且方差大致恒定。如果预测值越大、误差波动也越大,图中残差会出现扇形扩散,说明常方差假设被破坏了。第三个是样本之间的误差要相互独立,比如对同一家店铺重复抽样,误差可能不再独立。第四个是特征之间不要存在严重多重共线性,否则系数估计会很不稳定。

这些假设并不是用来吓唬人的,它们主要影响两件事:系数是否稳定、统计检验和置信区间是否可信。如果只是做简单预测,可能影响不大;但如果要拿系数给业务做解释,就必须回头检查这些前提。

2. 环境准备与模拟数据:先别急着上真实数据

2.1 本地运行需要的库与版本建议

开始写代码前,先确认环境。我平时用的组合是 numpy、pandas、matplotlib 和 scikit-learn,这四个库足够覆盖大部分回归场景。安装命令很简单:

pip install numpy pandas matplotlib scikit-learn

如果是在 Jupyter Notebook 里跑,我建议顺手把 numpy 和 sklearn 的版本打印出来,避免因为版本差异导致结果不一致。我个人在实操中见过不少因为 sklearn 版本更新、默认参数变化而导致结果对不上的案例,所以养成打印版本的习惯是值得的。

import numpy as np import pandas as pd import sklearn print(np.__version__) print(sklearn.__version__)

这一步看起来很小,但在复现别人的代码、或把代码分享给同事时,它能省掉非常多莫名其妙的沟通成本。

2.2 构造一份模拟数据,制造一个“已知答案”的测试场

学习回归时,我建议先不要直接上真实数据。真实数据里的噪声、缺失、共线性混在一起,出了问题你根本分不清是代码写错了,还是数据本身就长这样。先用模拟数据,等于先给自己建一个“已知答案”的测试场。

下面的代码会生成一份有三个变量的模拟数据集:两个特征 x1、x2,一个目标值 y。我故意把真实的系数和截距设置成已知数值,这样才能用来验证后面的代码算出来的结果是否接近真相。

import numpy as np import pandas as pd np.random.seed(42) n = 200 x1 = np.random.normal(0, 1, n) x2 = np.random.normal(0, 1, n) beta0_true = 3.2 beta1_true = 2.5 beta2_true = -1.8 noise_std = 1.0 y = beta0_true + beta1_true * x1 + beta2_true * x2 + np.random.normal(0, noise_std, n) df = pd.DataFrame({ "x1": x1, "x2": x2, "y": y })

为什么一定要加噪声?因为现实中几乎不存在“所有点正好落在一条直线上”的完美数据。加噪声也是在做一件符合直觉的事:模型只能解释一部分变化,剩下的是随机波动。这里的noise_std = 1.0控制了波动的幅度,你可以把它改大改小,感受一下噪声对拟合结果的影响。

模拟数据的另一个好处是,你可以随时把样本量 n 调大,比如改成 2000,然后观察估计出来的系数是不是越来越接近真实值。做一次这样的对照实验,比背十遍中心极限定理都管用。

2.3 数据生成过程带来的启发

构造完模拟数据后,先别急着进入下一步。可以把这份数据当作一面镜子,问自己几个问题:如果我现在用普通最小二乘法去拟合,系数能恢复到接近 2.5 和 -1.8 吗?如果我把噪声的标准差调到 3.0,系数估计值是不是明显变差?

实际在跑数据时,我也会想另一个问题:这份数据的生成过程到底是什么?你收集到销售数据,不代表销售就是曝光量、上架商品数两个变量决定的;很可能还有促销活动、天气、节假日等因素。当你建立的模型只包含部分变量时,实际上是在用一个简化模型去逼近一个复杂过程。理解这一点之后,你就不会对一个 R² 达到 0.95 的模拟数据过于兴奋,因为真实数据几乎不可能这么配合。

3. 手写最小二乘法:把回归原理落到数组运算上

3.1 正规方程的直观推导

要真正掌握线性回归,不能只会调用库函数,至少要能解释最小二乘解是怎么来的。我们将目标函数写成矩阵形式:

L(β) = (y - Xβ)^T (y - Xβ)

对 β 求导,并令导数为零:

∂L/∂β = -2X^T(y - Xβ) = 0

整理后得到:

X^T X β = X^T y

如果矩阵X^T X可逆,就有:

β = (X^T X)^(-1) X^T y

这套公式就是大家常说的正规方程。它最直观的优点是不需要像梯度下降那样手动选择学习率、不需要迭代,直接一步算出解析解。当特征维度不高、样本量不是特别大时,这是最简单粗暴的求解方式。

这里有一步很多人会忽略:不要把“线性回归”和“梯度下降”绑死。线性回归是模型,最小二乘是求参数的一种方法,梯度下降只是另一种求解方式。当数据量大到矩阵运算不划算时,梯度下降会更合适;但理解回归本身,还应该从正规方程这个起点开始。

3.2 用 numpy 实现一次最小二乘拟合

现在把公式落到代码里。我们先把包含截距项的设计矩阵构造出来,然后直接用 numpy 求解。

X = np.column_stack([ np.ones(n), # 截距项对应的常数特征 df["x1"].to_numpy(), df["x2"].to_numpy() ]) y = df["y"].to_numpy() # 正规方程:β = (XᵀX)⁻¹ Xᵀy beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y print("截距估计:", beta_hat[0]) print("x1系数估计:", beta_hat[1]) print("x2系数估计:", beta_hat[2])

如果一切正常,运行结果会非常接近我们预设的 3.2、2.5、-1.8。差异来自我们添加的随机噪声,只要噪声标准差不是太大,估计结果会在真值附近波动。

这里我强调一点:X.T @ X得到的是一个 3×3 的矩阵,np.linalg.inv求它的逆矩阵,再乘以X.T @ y,最终得到 β。写代码时可以一行完成,但建议初学者拆成两行,把中间矩阵和最终结果打印出来,理解每一步的维度变化。只有亲眼看到矩阵维度,才知道这一步背后的含义。

3.3 当正规方程不靠谱时

正规方程虽然直观,但在实际项目中往往不是首选求解方式。原因有几个:

第一个是计算复杂度。X^T X的求逆计算量大约是特征维度的三次方,当特征维度从几十变成几千、几万,正规方程会变得非常慢,甚至内存不够。第二个是数值稳定性。如果特征之间存在高度相关性,X^T X会接近奇异,求逆容易得到极大或极不稳定的结果。第三个是面对稀疏数据时,正规方程会把原本可以用稀疏结构解决的问题变得笨重。

所以,更合适的做法是使用np.linalg.lstsq,它在底层用了更稳定的数值算法。对同一份数据,调用方式如下:

beta_hat_stable, residuals, rank, s = np.linalg.lstsq(X, y, rcond=None) print("更稳定的系数估计:", beta_hat_stable)

如果你手动写回归代码只是为了学习,我建议用正规方程跑通一次,然后记住一句话:真正生产代码里,除非特征数很少且矩阵条件良好,否则优先用lstsq或封装好的库。

4. 使用 sklearn:调用工具前的几个关键选择

4.1 最小可用代码与训练测试拆分

手写版本跑通之后,就可以放心用封装工具提升效率了。sklearn 的LinearRegression接口很简单,但有几处细节需要留意。先看完整流程:

from sklearn.model_selection import train_test_split from sklearn.linear_model import LinearRegression X = df[["x1", "x2"]] y = df["y"] X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.2, random_state=42 ) model = LinearRegression() model.fit(X_train, y_train) print("截距:", model.intercept_) print("系数:", model.coef_) y_pred = model.predict(X_test)

很多第一次接触 sklearn 的人会把重心放在fit上,但train_test_split这一步同样重要。如果直接把所有数据拿去fit,你无法知道模型对新数据的泛化能力。保留一部分测试数据,相当于考试时先把答案遮住一部分,最后再拿出来验证。

在LinearRegression的默认参数里,fit_intercept=True。这意味着 sklearn 会自动帮你处理截距项,不需要像手写时那样手动构造一列 1。这个默认行为很方便,但如果你之前已经手动给数据加过一列常数特征,再设置fit_intercept=True,就会出现完全共线性,导致结果出问题。可以记住:二选一,不要重复加截距。

4.2 归一化:什么时候需要,什么时候不需要

关于线性回归是否要做归一化,网上说法很多。我直接说结论:普通最小二乘对特征是否归一化其实不敏感,因为系数会等比例缩放来抵消量纲差异;真正需要注意的是正则化回归,比如岭回归和 Lasso。

那为什么很多教程强调要归一化?主要是为了可解释性和正则化的公平性。当两个特征量纲差异极大,比如一个特征是“曝光量”以万为单位,另一个特征是“转化率”以百分比为单位,你得到的原始系数大小不可直接比较。归一化之后,每个特征都变成均值为 0、标准差为 1 的量,系数就能从“每增加一个单位”变成“每增加一个标准差”。

如果你打算在回归里加 L1 或 L2 惩罚,那归一化几乎是必须的。因为惩罚项对所有系数一视同仁,如果特征量纲不同,惩罚就会被量纲大的特征主导。

用 sklearn 的StandardScaler配合管道,可以避免很多手动作业:

from sklearn.preprocessing import StandardScaler from sklearn.pipeline import make_pipeline pipeline = make_pipeline(StandardScaler(), LinearRegression()) pipeline.fit(X_train, y_train) y_pred_scale = pipeline.predict(X_test) print("缩放后的系数:", pipeline.named_steps["linearregression"].coef_) print("缩放后的截距:", pipeline.named_steps["linearregression"].intercept_)

使用管道时,StandardScaler会自动用训练集的数据计算均值和标准差,再对训练集和测试集分别做相同的变换,不会把测试集信息泄漏到训练过程里。

4.3 系数解读:从数字回到业务语言

模型训练完之后,不要只看 R²,还要能解释系数。在小店营业额的例子里,如果 x1 是广告曝光量,它的系数是 2.5,意思是“在其他条件不变的情况下,广告曝光量每增加一个单位,预测营业额平均增加 2.5 个单位”。

注意“其他条件不变”这个前提。真实数据中特征之间总是相互关联,很难真正控制其他变量不变,但这是线性回归给出可解释方向的基本语境。如果特征经过了标准化,这个含义要改成“广告曝光量每增加一个标准差,预测营业额平均增加对应的多少个单位”。

系数解释中最容易翻车的是符号方向。如果模型给出的系数符号和业务直觉相反,不能强行解释,需要回到第 6 部分排查原因,大概率是共线性或数据泄露,而不是业务逻辑错了。

5. 评价模型:看 R² 远远不够,残差才说实话

5.1 几个常用指标的计算和解读

训练完后需要用指标评估。回归问题最常用的三个指标是均方误差(MSE)、平均绝对误差(MAE)和决定系数(R²)。用 sklearn 可以一次算出来:

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score mse = mean_squared_error(y_test, y_pred) mae = mean_absolute_error(y_test, y_pred) r2 = r2_score(y_test, y_pred) print("MSE:", mse) print("MAE:", mae) print("R²:", r2)
指标公式含义特点适合场景
MSE预测误差的平方取平均大误差会被放大希望避免严重偏差时
MAE预测误差的绝对值取平均受异常值影响较小对离群点容忍度更高时
R²1 减去残差平方和占总平方和的比例越接近 1 越好,可能为负看模型的相对解释能力

R² 很容易让人产生误解。模拟数据里 R² 可以到 0.95,但真实业务数据 R² 等于 0.3 也很常见,这不代表模型没有用。R² 只说明模型解释了目标变量多少比例的变化,它本身不能判断预测误差是否能被业务接受。

另一个容易踩的坑是 MSE 的量纲。它跟目标变量的平方是同一个量纲,所以 1000 的 MSE 到底大不大,取决于营业额本身是几十还是几万。A 任务的 MSE 是 100,可能已经很好;B 任务的 MSE 是 10,可能仍然不可接受。

5.2 残差分析:画一张散点图往往比卡指标更有用

指标是把所有误差压缩成一个数,而残差图能让你看到误差的结构。最简单的做法是画“预测值 vs 残差”的散点图:

import matplotlib.pyplot as plt residual = y_test - y_pred plt.scatter(y_pred, residual, alpha=0.6) plt.axhline(y=0, color="gray", linestyle="--") plt.xlabel("预测值") plt.ylabel("残差") plt.show()

如果模型拟合得当,残差应该大致随机分布在水平线 0 的两侧,没有明显的漏斗形或曲线形。如果看到残差随着预测值增大而逐渐发散,说明方差不是常数,最常用的对策是对目标变量取对数后再建模。如果残差整体呈现抛物线形状,说明某个特征和目标的真实关系是非线性的,模型中缺少一个平方项或交互项。

残差图看起来像一个探索过程,而不是一个“看结果好不好”的关卡。我会在真实项目里先画图,再决定要不要做特征变换,而不是上来就把所有指标调一遍。

5.3 训练集与测试集的配合

对同一份数据,应该同时看训练集和测试集的结果。一个常见的现象是训练 R² 远高于测试 R²,说明模型把训练数据背下来了,也就是过拟合。

线性回归在特征很少时通常不容易过拟合,但一旦特征数量接近样本数量,或存在强共线性,过拟合就不可避免。为了避免只靠一次随机划分下结论,可以在代码里加入交叉验证:

from sklearn.model_selection import cross_val_score scores = cross_val_score(model, X, y, cv=5, scoring="neg_mean_squared_error") print("负MSE均值:", scores.mean())

这里 scikit-learn 用neg_mean_squared_error作为评分,因为cross_val_score默认认为分数越大越好,所以对误差取负数。看到负数不要慌,这只是评分方向问题,实际想看 MSE 时取一次绝对值即可。

6. 实战中最常见的几个坑和排查思路

6.1 系数符号与业务常识完全相反

这是我在实际项目里遇到最多的问题,没有之一。小例子是这样的:理论上投入越多产出应该更高,但模型给出的系数是负的,而且绝对值还很大。

出现这个现象的原因通常有三个。第一是多重共线性,x1 和 x2 高度相关,模型无法稳定地拆解出每个特征的独立贡献,某个系数就会被另一个特征“抢走”一部分解释力。第二是漏掉了关键变量,导致某个重要因素的作用被错误地归到现有特征头上。第三是数据本身就存在反向因果关系,比如营业额高的月份广告曝光也高,但真正驱动营业额的是季节性,而不是广告投放本身。

排查思路不要上来就改特征,而是先做最小实例验证:只把 x1 和目标变量做一次单变量回归,看系数方向。如果单变量方向正常,加入其他变量后方向变了,多半是共线性或混杂。再打印特征之间的相关矩阵,看看有没有 |r| 接近 0.9 或更高的特征对。

6.2 多重共线性导致系数不稳定

共线性真正麻烦的地方在于,模型整体预测效果可能还不错,但单个系数变得极其不稳定。换一两行数据,系数就可能从 2 变成 -1。这种模型没法跟业务解释。

判断多重共线性,最直接的方式是看特征相关矩阵。如果两个特征的皮尔逊相关系数绝对值大于 0.8,就要引起警惕。更严格的做法是计算方差膨胀因子,但实际操作中我会先做相关矩阵,因为它更直观。

处理方式大致有三种:一是直接删除高度相关的特征之一,简单有效,适合特征本来就不多的情况;二是用岭回归或 Lasso 这类带正则化的模型,让系数收缩,换来稳定性;三是把相关特征做 PCA 降维后再建模,但这会让系数解释变得复杂。大部分业务场景我会优先用“删除或合并特征”的方案。

6.3 数据泄露让模型成绩虚高

数据泄露是另一个绕不开的坑。它的本质是,模型在训练时接触到了测试阶段不该出现的信息,导致测试成绩失真。最常见的泄露发生在预处理阶段。

比如有人先对全量数据做标准化,再拆分训练集和测试集。标准化用到的均值和标准差来自全量数据,其中包含了测试集的信息。虽然这种泄露对线性回归的影响通常不如对神经网络大,但它仍然会让评估结果偏乐观。

# 错误示范:先缩放再拆分 from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_scaled = scaler.fit_transform(X) X_train, X_test, y_train, y_test = train_test_split(X_scaled, y)

正确做法是先拆分,再对训练集执行fit_transform,对测试集只用transform,或者更省事地直接放进Pipeline,让 sklearn 自己控制顺序。

更隐蔽的泄露是特征构造。比如预测“下个月销售额”时,特征里不小心混入了“本月的实际销售额”,测试集上模型成绩会非常漂亮,但到了真正要预测未来时彻底失效。

6.4 我的问题排查顺序

现在整理一下我遇到的回归问题时的排查顺序,这个顺序适用性很广:

  • 第一步想象数据本身:哪些特征有缺失?哪些特征取值范围明显异常?是否用错了单位?
  • 第二步检查拆分:数据在 fit 之前是否已经被全量统计量污染?
  • 第三步看特征关系:打印相关矩阵,找出高度相关特征,做一轮单变量回归。
  • 第四步看残差图:确认残差是否随机分布,判断是否需要做非线性变换。
  • 第五步对比训练集和测试集分数:差距非常大时优先怀疑过拟合和数据泄露。

这个顺序不一定解决所有问题,但至少能把“代码 bug”“数据问题”“模型方法问题”三类原因分开,避免在某个环节上反复打转。

7. 下一步升级:从普通线性回归走向更复杂的选择

7.1 正则化回归:岭回归和 Lasso 的直觉

当特征变多、共线性明显时,普通线性回归的系数会越来越大,模型变得脆弱。正则化回归的思路是在原来的损失函数后面追加一个惩罚项,逼着系数不要疯涨。

岭回归加的是 L2 惩罚,对系数做平方和约束,倾向于把不太重要的系数压到接近零,但不会精确等于零。Lasso 加的是 L1 惩罚,对系数绝对值求和,能把一部分系数精确变成零,相当于顺手做了特征选择。两者的使用场景略有侧重,但代码写法非常接近:

from sklearn.linear_model import Ridge, Lasso ridge = Ridge(alpha=0.5) ridge.fit(X_train, y_train) lasso = Lasso(alpha=0.01) lasso.fit(X_train, y_train) print("Ridge系数:", ridge.coef_) print("Lasso系数:", lasso.coef_)

这里的 alpha 是正则化强度,越大惩罚越重。实际选 alpha 时不要手调,用RidgeCV、LassoCV这类带交叉验证的版本自动搜索更好。

7.2 特征变换:当线性关系不够时

现实数据里,目标变量和特征的真实关系很少是完美的直线。比如“面积”对“房价”的影响,可能是边际递减的;广告曝光量对销售额的促进作用,也可能存在饱和效应。

最简单的手段是给特征加多项式项和交互项。sklearn 里有一个很好用的类:

from sklearn.preprocessing import PolynomialFeatures poly = PolynomialFeatures(degree=2, include_bias=False) X_poly = poly.fit_transform(X) model_poly = LinearRegression() model_poly.fit(X_poly, y)

degree=2 会把原来的每个特征扩展成平方项,并自动生成两两交互项。在数据量不大、特征不多时,这是快速提升拟合能力的手段。但也要注意,多项式特征会让特征数量快速膨胀,过拟合风险随之上升,这时可以搭配岭回归或 Lasso 使用,效果通常比单独用更靠谱。

7.3 我自己的实操体会

写到这里,我不太想按惯例做一个“总结”,反而想留下一个最实际的经验。

我现在的习惯是:无论面对什么回归任务,都会先手写最小二乘跑一遍,再对比 sklearn 的结果。这个习惯表面上是在做双重验证,实际上是在把思考锚定在一个非常具体的点上:我的数据和特征长什么样?样本量是否够?特征之间是否独立?如果你能对这组问题给出清晰答案,那用哪个库都只是工具差异,你的工作重心会从“调库调参”变成“理解问题”。

别把线性回归当成入门后就要扔掉的玩具。很多复杂模型的输出本质上包含了线性结构;异常值检测、因果推断里的倾向得分、时间序列里的自回归,处处都有它的影子。把这一块地基打牢,后面学任何模型都会轻松很多。

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进