最小二乘法：从误差量化到多元回归，原理推导与Python实践

1. 从“猜”到“算”：为什么我们需要最小二乘法？

做数据分析、搞模型拟合，或者哪怕只是用Excel画条趋势线，你大概率都听过“最小二乘法”这个名字。它听起来像个高深的数学工具，但实际上，它的核心思想朴素得惊人：找一条线（或者一个面），让所有数据点到这条线的“距离”平方和最小。

想象一个场景：你是个质量工程师，要研究生产线上某个零件的加工时间（X）和最终尺寸误差（Y）的关系。你测了10组数据，点在坐标图上散落一片。老板问你：“这俩变量到底啥关系？能不能用个公式大致描述一下？”你当然可以凭感觉画条直线穿过去，但张三画的和你画的可能不一样，谁对谁错？缺乏一个客观标准。

或者，你是个开发者，在做一个智能家居的温控模型。传感器传回的温度数据和设定值之间总有偏差，你需要一个算法来自动调整加热功率，使得实际温度尽可能平稳地跟随设定曲线。这个“尽可能”怎么量化？怎么让机器自动找到最优的调整参数？

最小二乘法就是解决这类问题的“标准答案”和“自动化工具”。它不靠猜，而是靠算。通过严谨的数学推导，它给出了一组公式，只要你把数据代进去，就能唯一地、最优地确定那条拟合曲线的参数。这个“最优”的标准，就是前面说的“距离的平方和最小”。这个思想在两百多年前由高斯和勒让德等人奠定，至今仍是科学计算、工程分析、机器学习（尤其是线性回归）的基石。

我最初接触它时，觉得那一堆求和符号和偏导数很吓人。但后来在调参、分析实验数据、甚至理解一些机器学习算法内部运作时，一次次重温，才发现它的精妙之处不在于复杂的推导，而在于它用非常简洁的数学目标（最小化平方和），优雅地解决了无数实际中的近似和估计问题。接下来，我会抛开那些让人望而生畏的教科书式证明，带你从问题出发，一步步拆解最小二乘法的核心逻辑、手算与代码实现、以及那些容易踩坑的细节，让你不仅能用它，更能懂它为什么这么用。

2. 核心思想拆解：误差、平方和与“最小”的奥义

要理解最小二乘法，关键在于吃透它的三个核心概念：误差定义、为什么用平方、以及“最小化”如何实现。我们以最简单也最常用的一元线性回归为例，即用一条直线y = ax + b来拟合数据。

2.1 如何定义“不准”？—— 误差的量化

假设我们有n个数据点(x_i, y_i),i = 1, 2, ..., n。我们用一条直线ŷ_i = a * x_i + b去预测每一个x_i对应的y值。这里ŷ_i（读作 y-hat）表示预测值。

那么，对于第i个点，预测值ŷ_i和真实值y_i之间的差距，就是误差（Error）或残差（Residual），记作e_i：e_i = y_i - ŷ_i = y_i - (a * x_i + b)

这个e_i可正可负。如果点在线之上，误差为正；点在线之下，误差为负。

2.2 为什么是“平方”和？—— 克服正负相消与放大大误差

现在我们有了每个点的误差。如何评价整条直线y = ax + b的好坏？一个朴素的想法是把所有误差加起来：S = e_1 + e_2 + ... + e_n。但这会出大问题：正误差和负误差会相互抵消！一条严重偏离的直线，如果其误差正负平衡，总和可能接近零，这显然不合理。

于是，我们想到用绝对值：S_abs = |e_1| + |e_2| + ... + |e_n|。这确实避免了抵消问题，从数学上讲，最小化绝对值和（L1范数）也是完全可行的，并且有它的优势（如对异常值更鲁棒）。但在历史上和大多数实际应用中，最小二乘法选择了平方和。

选择平方和的几个关键理由：

1. 数学处理友好：绝对值函数在零点不可导，而平方函数处处光滑可导。这使得我们可以使用强大的微积分工具（求导）来寻找最小值点，过程会简洁优雅得多。
2. 强调大误差：平方操作会放大较大的误差。例如，一个误差为2的点贡献为4，而一个误差为10的点贡献高达100。这意味着最小二乘法会极力避免出现大的偏差，拟合出的直线会尽可能靠近那些误差较大的点，从而在整体上更均衡。这在许多工程和科学场景中是 desired 的特性。
3. 统计意义：在高斯-马尔可夫定理的假设下（误差零均值、同方差、不相关），最小二乘估计是所有线性无偏估计中方差最小的（BLUE, Best Linear Unbiased Estimator）。这为它提供了坚实的统计学基础。

因此，我们定义目标函数（损失函数）Q为所有误差的平方和：Q(a, b) = Σ(e_i)^2 = Σ [y_i - (a * x_i + b)]^2， 其中Σ表示对i从1到n求和。

我们的目标变得非常清晰：找到一对参数(a, b)，使得这个平方和Q(a, b)达到最小。

2.3 如何找到“最小”点？—— 微积分的威力

现在问题转化成了一个二元函数Q(a, b)的优化问题。回想微积分知识：对于一个光滑函数，在其极小值点处，它对各个自变量的偏导数应为零。

所以我们分别对a和b求偏导，并令其等于零：

1. 对b求偏导：∂Q/∂b = Σ 2 * [y_i - (a*x_i + b)] * (-1) = -2 Σ [y_i - a*x_i - b] = 0化简得：Σ y_i - a Σ x_i - n * b = 0-->n*b + a Σ x_i = Σ y_i。(方程1)

2. 对a求偏导：∂Q/∂a = Σ 2 * [y_i - (a*x_i + b)] * (-x_i) = -2 Σ [x_i * (y_i - a*x_i - b)] = 0化简得：Σ (x_i * y_i) - a Σ (x_i^2) - b Σ x_i = 0-->b Σ x_i + a Σ (x_i^2) = Σ (x_i * y_i)。(方程2)

这样我们得到了关于未知数a和b的正规方程组（Normal Equations）：

n*b + (Σ x_i) * a = Σ y_i (Σ x_i) * b + (Σ x_i^2) * a = Σ (x_i * y_i)

这是一个二元一次方程组。解这个方程组，就能得到a和b的解析解（公式解）：

a = [n * Σ(x_i*y_i) - Σ x_i * Σ y_i] / [n * Σ(x_i^2) - (Σ x_i)^2]

b = [Σ y_i - a * Σ x_i] / n = ȳ - a * x̄

其中，x̄和ȳ分别是x和y的样本均值。b的公式非常直观：最优直线必然穿过数据的中心点(x̄, ȳ)。

注意：这里有一个经典的“坑”。计算a的分母n * Σ(x_i^2) - (Σ x_i)^2，在数学上它等于n * Σ (x_i - x̄)^2，即x的方差乘以(n-1)。这个分母必须不为零。如果为零，意味着所有x_i都相同，即数据点在一条竖直线上，此时x的方差为零，不存在唯一的斜率a。这在物理意义上也很好理解：当x值没有变化时，我们无法衡量y随x的变化率。

至此，我们从定义问题（量化误差），到选择优化目标（最小化平方和），再到运用数学工具（求导解方程），完整地推导出了一元线性最小二乘法的核心公式。这个过程本身，就是理解其精髓的关键。

3. 从公式到代码：手算验证与Python实现

理解了原理，我们最好通过一个具体的例子来巩固。手动计算能加深对公式每个部分的理解，而代码实现则是将其应用于实际问题的必经之路。

3.1 手动计算示例

假设我们有以下5组数据，研究学习时间(X)与考试成绩(Y)的关系：

学习时间 (x)	考试成绩 (y)	x*y	x^2
2	65	130	4
4	75	300	16
6	85	510	36
8	95	760	64
10	105	1050	100
Σ	Σx=30	Σy=425	Σxy=2750

这里 n=5。我们代入公式计算：

1. 计算斜率a： 分子 =n * Σ(xy) - Σx * Σy = 5*2750 - 30*425 = 13750 - 12750 = 1000分母 =n * Σ(x^2) - (Σx)^2 = 5*220 - 30^2 = 1100 - 900 = 200所以a = 1000 / 200 = 5

2. 计算截距b： 先求均值：x̄ = 30/5 = 6,ȳ = 425/5 = 85b = ȳ - a * x̄ = 85 - 5*6 = 85 - 30 = 55

因此，我们得到拟合直线为：ŷ = 5*x + 55。

解读：斜率a=5意味着，平均来说，学习时间每增加1小时，考试成绩预计提高5分。截距b=55可以理解为，当学习时间为0时，基础的预期成绩（可能包含了其他因素或基础能力）。你可以将x值代回方程，计算预测值ŷ，会发现它们完美地落在一条直线上（因为本例数据是人为构造的线性关系）。

3.2 Python代码实现：三种常用方法

在实际工作中，我们几乎不会手动计算，而是借助工具。Python的NumPy和SciPy库提供了极其便捷的接口。下面展示三种主流方法。

方法一：使用NumPy的polyfit函数（最简洁）polyfit是专门用于多项式拟合的函数，一元线性拟合就是一次多项式拟合。

import numpy as np # 数据 x = np.array([2, 4, 6, 8, 10]) y = np.array([65, 75, 85, 95, 105]) # 进行1次多项式（线性）拟合，返回系数，最高次幂在前 coefficients = np.polyfit(x, y, deg=1) # coefficients[0]是斜率a, coefficients[1]是截距b a_np, b_np = coefficients print(f"NumPy polyfit 结果: 斜率 a = {a_np:.2f}, 截距 b = {b_np:.2f}") # 输出: 斜率 a = 5.00, 截距 b = 55.00

方法二：使用NumPy的线性代数方法（理解正规方程）这种方法直接求解我们之前推导的正规方程组(X^T * X) * β = X^T * y，其中X是设计矩阵（第一列为1，第二列为x），β是参数向量[b, a]^T。

import numpy as np x = np.array([2, 4, 6, 8, 10]) y = np.array([65, 75, 85, 95, 105]) # 构造设计矩阵 X: 第一列全1（对应截距b），第二列为x值 X = np.vstack([np.ones_like(x), x]).T # .T表示转置 # 现在 X 是一个 5行2列的矩阵，每行是 [1, x_i] # 求解正规方程: β = (X^T * X)^(-1) * X^T * y # 使用 np.linalg.inv 求逆，@ 表示矩阵乘法 beta = np.linalg.inv(X.T @ X) @ X.T @ y b_la, a_la = beta # 注意顺序，beta[0]是截距，beta[1]是斜率 print(f"线性代数方法 结果: 截距 b = {b_la:.2f}, 斜率 a = {a_la:.2f}") # 输出: 截距 b = 55.00, 斜率 a = 5.00

方法三：使用SciPy的stats.linregress（功能更全）scipy.stats.linregress不仅返回参数，还提供丰富的统计量，如R值（相关系数）、p值、标准误等，非常适合统计分析。

from scipy import stats x = [2, 4, 6, 8, 10] y = [65, 75, 85, 95, 105] # 执行线性回归 slope, intercept, r_value, p_value, std_err = stats.linregress(x, y) print(f"SciPy linregress 结果:") print(f" 斜率 a = {slope:.2f}") print(f" 截距 b = {intercept:.2f}") print(f" 相关系数 R = {r_value:.4f}") print(f" R^2 = {r_value**2:.4f}") # 决定系数，拟合优度 print(f" 斜率标准误 = {std_err:.4f}") print(f" p值 = {p_value:.4f}") # 输出中 R=1.0, p值极小，说明拟合极好。

实操心得：

1. 日常快速拟合，用np.polyfit最方便。
2. 需要深入理解矩阵运算或自定义更复杂的模型（如多元、带权重），方法二（正规方程）的思路是基础。但注意，对于特征非常多（列数很大）的情况，直接求逆np.linalg.inv(X.T @ X)可能数值不稳定或计算慢，此时会采用梯度下降等迭代法。
3. 需要进行统计推断（看显著性、置信区间等），scipy.stats.linregress或更专业的statsmodels库是更好的选择。
4. 一个常见坑点：数据量纲。如果x是“万像素”，y是“亿元销售额”，计算出的斜率a会非常大且难以解释。通常建议对数据进行标准化（减均值除以标准差）处理，这样得到的斜率是“变化一个标准差x，引起y变化多少个标准差”，更具可比性。

4. 多元线性回归：当影响因素不止一个

现实问题中，结果变量y往往受多个因素x1, x2, ..., xp影响。例如，房价可能受面积、卧室数量、房龄、地段等多个因素影响。这时，我们就需要将一元情况推广到多元线性回归，模型变为：ŷ = b0 + b1*x1 + b2*x2 + ... + bp*xp其中b0是截距，b1到bp是各自变量对应的系数。

最小二乘法的思想完全不变：寻找一组参数(b0, b1, ..., bp)，使得预测值ŷ与真实值y之差的平方和最小。

4.1 矩阵形式：优雅的统一

使用矩阵表示会异常简洁。令：

• y是一个n×1的列向量，包含所有观测值[y1, y2, ..., yn]^T。
• X是一个n×(p+1)的设计矩阵。它的第一列全是1（对应截距b0），后面p列分别是各个自变量x1, x2, ..., xp的观测值。
• β是一个(p+1)×1的列向量，包含所有待求参数[b0, b1, ..., bp]^T。
• ε是一个n×1的列向量，表示误差。

那么，整个模型可以写成：y = Xβ + ε

最小二乘的目标函数Q(β) = Σ(y_i - ŷ_i)^2可以写成向量形式：Q(β) = (y - Xβ)^T (y - Xβ)

通过对β求导（向量求导），并令导数为零向量，我们可以得到正规方程组的矩阵形式：X^T X β = X^T y

如果X^T X是可逆矩阵（要求X列满秩，即自变量之间不存在严格的线性相关），那么最优参数向量的解析解为：β = (X^T X)^{-1} X^T y

这个公式在形式上和一元的解a = ...惊人地统一，体现了矩阵表达的威力。在Python中，我们依然可以用方法二（np.linalg.inv）来求解，但更稳健的做法是使用np.linalg.lstsq或np.linalg.solve。

4.2 Python实现与解读

假设我们研究房价(price)，考虑面积(area)和卧室数(bedrooms)两个因素。数据如下：

area	bedrooms	price
100	2	300
150	3	450
120	2	320
180	4	550
90	1	260

import numpy as np # 数据 X_features = np.array([[100, 2], [150, 3], [120, 2], [180, 4], [90, 1]]) # 特征矩阵，不包含截距列 y = np.array([300, 450, 320, 550, 260]) # 方法：使用 np.linalg.lstsq 求解最小二乘解（推荐，数值更稳定） # 首先，需要为X添加一列1，用于估计截距 X_design = np.c_[np.ones(X_features.shape[0]), X_features] # 在左侧添加一列1 # 使用 lstsq 求解，它通过SVD分解求解，比直接求逆稳定 beta, residuals, rank, s = np.linalg.lstsq(X_design, y, rcond=None) # beta 包含了 b0, b1, b2 b0, b1, b2 = beta print(f"多元回归结果:") print(f" 截距 b0 = {b0:.2f}") print(f" 面积系数 b1 = {b1:.2f}") print(f" 卧室系数 b2 = {b2:.2f}") print(f" 模型公式: price = {b0:.2f} + {b1:.2f}*area + {b2:.2f}*bedrooms") # 计算预测值 y_pred = X_design @ beta print(f" 预测房价: {y_pred}")

运行后，你可能会得到类似price = 50.12 + 2.05*area + 25.88*bedrooms的结果。

系数解读：

• b1=2.05：在卧室数量保持不变的情况下，面积每增加1平方米，房价平均上涨约2.05（单位）。
• b2=25.88：在面积保持不变的情况下，卧室每增加1间，房价平均上涨约25.88（单位）。
• b0=50.12：当面积和卧室数都为0时的基础房价（这个解释在现实中可能无实际意义，更多是数学上的截距）。

重要注意事项（多元回归的坑）：

1. 多重共线性：这是多元回归中最常见也最棘手的问题之一。如果两个自变量（如“面积”和“卧室数”）高度相关，那么X^T X矩阵会接近奇异（不可逆），导致系数估计(X^T X)^{-1}变得极不稳定，系数方差很大，解释性变差。表现为：系数符号不符合常识、微小的数据变动导致系数巨大变化。诊断方法：计算方差膨胀因子(VIF)。解决方法：剔除相关性高的变量之一、使用主成分回归(PCR)或岭回归(Ridge Regression)等正则化方法。
2. 特征缩放：当自变量的量纲和数值范围差异巨大时（如“面积”是100-200，“家庭年收入”是10万-50万），未经缩放的梯度下降算法可能收敛很慢，且正则化惩罚项会对大数值特征产生不公平的影响。虽然对于解析解(X^T X)^{-1} X^T y来说，缩放不影响最终预测，但会影响系数的解释。通常建议进行标准化（Standardization）或归一化（Normalization）。
3. 过拟合：当自变量数量p很多，甚至接近样本量n时，模型很容易完美拟合训练数据（平方和接近零），但在新数据上表现很差。这需要通过训练集-测试集划分、交叉验证来评估，并使用正则化（如Lasso, Ridge）或特征选择来缓解。

5. 非线性关系的处理：多项式回归与线性化

最小二乘法本质上是“线性”的，这个“线性”指的是参数是线性的，即模型关于待求参数β是线性的。但并不意味着y和x的关系必须是直线。我们可以通过巧妙的变换，将许多非线性关系转化为线性模型来处理。

5.1 多项式回归

如果y和x的关系是曲线，例如二次、三次关系，我们可以令x1 = x,x2 = x^2,x3 = x^3，然后拟合模型：ŷ = b0 + b1*x + b2*x^2 + b3*x^3这依然是一个关于参数b0, b1, b2, b3的线性模型！我们只需要将原始特征x扩展为多项式特征[x, x^2, x^3]，然后套用多元线性回归的框架即可。

Python实现（使用sklearn）：

import numpy as np from sklearn.preprocessing import PolynomialFeatures from sklearn.linear_model import LinearRegression import matplotlib.pyplot as plt # 生成非线性数据 np.random.seed(42) x = np.linspace(-3, 3, 100) y_true = 0.5 * x**2 + x + 2 # 真实的二次关系 y_noise = y_true + np.random.randn(100) * 1.5 # 加入噪声 # 将数据转换为二维数组（sklearn要求） x_reshaped = x.reshape(-1, 1) # 创建多项式特征（最高2次） poly = PolynomialFeatures(degree=2, include_bias=False) # include_bias=False，因为LinearRegression自带截距 X_poly = poly.fit_transform(x_reshaped) # 现在X_poly有两列：[x, x^2] # 使用线性回归拟合 model = LinearRegression() model.fit(X_poly, y_noise) # 查看系数 print(f"截距: {model.intercept_:.4f}") print(f"系数 (对应 [x, x^2]): {model.coef_}") # 预测并绘图 y_pred = model.predict(X_poly) plt.scatter(x, y_noise, s=10, alpha=0.6, label='原始数据（含噪声）') plt.plot(x, y_true, 'r-', label='真实关系 (y=0.5x^2+x+2)') plt.plot(x, y_pred, 'g--', linewidth=2, label='二次多项式拟合') plt.legend() plt.xlabel('x') plt.ylabel('y') plt.title('多项式回归示例') plt.show()

你会发现，拟合出的系数[b1, b2]接近[1, 0.5]，截距接近2，成功还原了真实的数据生成过程。

注意：多项式阶数degree不宜过高。过高的阶数会导致模型过于复杂，在数据点之间剧烈震荡，即过拟合。需要通过交叉验证来选择合适的多项式次数。

5.2 可线性化的非线性模型

有些经典的非线性模型，可以通过简单的变量代换，转化为线性形式。

1. 指数模型：y = a * e^(b*x)

   • 线性化方法：两边取自然对数，ln(y) = ln(a) + b*x。
   • 令Y' = ln(y),A = ln(a)，则模型变为Y' = A + b*x，对(x, ln(y))做线性拟合即可。

2. 幂律模型：y = a * x^b

   • 线性化方法：两边取对数（常用以10为底或自然对数），log(y) = log(a) + b * log(x)。
   • 令Y' = log(y),X' = log(x),A = log(a)，则模型变为Y' = A + b*X'，对(log(x), log(y))做线性拟合。

3. 对数模型：y = a + b * ln(x)

   • 线性化方法：直接令X' = ln(x)，模型变为y = a + b*X'，对(ln(x), y)做线性拟合。

重要提醒：在对y进行变换（如取对数）后，我们最小化的是变换后变量的误差平方和（例如Σ(ln(y_i) - ln(ŷ_i))^2），而不是原始y的误差平方和。这会导致拟合目标不同。通常，如果误差结构在变换后的尺度上更接近正态分布、方差更稳定，这种方法是合适的。否则，更好的方法是使用非线性最小二乘法，直接最小化原始y的误差平方和，但这需要迭代优化算法（如梯度下降、Levenberg-Marquardt）。

6. 评估与诊断：你的模型真的靠谱吗？

拟合出参数只是第一步，评估模型好坏至关重要。不能只看“拟合线画上去挺好看”。

6.1 核心评估指标

1. 残差图（Residual Plot）：这是最直观、最重要的诊断工具。绘制预测值ŷ（或自变量x）与残差e_i = y_i - ŷ_i的散点图。

   • 理想情况：残差随机、均匀地分布在0线上下，没有明显的规律或趋势（如下图左）。
   • 出现问题：
     • 漏斗形：残差随ŷ增大而散开，提示异方差性，误差方差不是常数。
     • 曲线形：残差呈现U型或倒U型分布，提示模型非线性关系未被捕捉，可能需要加入高次项或交互项。
     • 离群点：个别点残差绝对值远大于其他点，可能是异常值，需要检查。

2. 决定系数 R²：最常用的拟合优度指标。R² = 1 - (SS_res / SS_tot)。

   • SS_res是残差平方和Σ(y_i - ŷ_i)^2。
   • SS_tot是总平方和Σ(y_i - ȳ)^2，反映了y自身的波动。
   • R²表示模型能够解释的y的方差比例，范围在0到1之间。越接近1，说明模型解释力越强。
   • 注意：R²会随着自变量增加而自然增大，即使加入无关变量。因此对于多元回归，更常用调整后R²，它惩罚了变量个数。

3. 均方误差（MSE）与均方根误差（RMSE）：

   • MSE = SS_res / n。数值越小越好，但它的大小依赖于y的量纲。
   • RMSE = sqrt(MSE)。它与y同量纲，更易于解释。例如，房价预测的RMSE是5万元，可以直观理解为“平均预测误差在5万元左右”。

4. 系数显著性检验（t检验）：在统计框架下，我们关心每个自变量x_j的系数b_j是否显著不为零。这通过计算t = b_j / SE(b_j)（其中SE是标准误）并查t分布表得到p值。通常p < 0.05认为该变量对y有显著影响。scipy.stats.linregress和statsmodels会提供这些统计量。

6.2 常见问题与对策

• 异方差性：残差方差不等。这违背了最小二乘法的经典假设，会导致系数标准误估计不准确，进而影响显著性检验。对策：加权最小二乘法(WLS)、对因变量进行变换（如取对数）、使用稳健标准误。
• 自相关性：时间序列数据中，残差前后相关。这也会导致标准误低估。对策：检查Durbin-Watson统计量，使用时间序列模型（如ARIMA）或广义最小二乘法(GLS)。
• 异常值与高杠杆点：异常值（y值异常）和高杠杆点（x值异常）会严重扭曲拟合结果。对策：绘制库克距离（Cook‘s Distance）图，识别对模型影响过大的点，并检查其数据是否正确，或考虑使用稳健回归方法（如RANSAC, Theil-Sen）。
• 模型设定错误：遗漏重要变量、或包含了无关变量、或函数形式错误（该用非线性却用了线性）。对策：基于领域知识选择变量，利用残差图、偏回归图诊断，尝试不同的模型形式。

一个简单的诊断流程：

1. 拟合模型后，首先绘制残差图，观察是否随机分布。
2. 查看R²和RMSE，对模型整体解释力和预测误差有个大致判断。
3. 对于多元回归，查看系数估计值、标准误、t值和p值，判断每个变量的显著性和影响方向。
4. 检查方差膨胀因子(VIF)，诊断多重共线性（通常VIF>10认为存在严重共线性）。
5. 计算库克距离，检查是否有强影响点。

7. 超越普通最小二乘法：正则化与稳健回归

当数据存在我们上面提到的一些问题时，普通最小二乘法（OLS）可能不是最优选择。现代统计学和机器学习提供了许多增强版本。

7.1 正则化：岭回归与Lasso回归

当自变量很多、存在多重共线性或为了防止过拟合时，我们可以在损失函数中加入一个对系数大小的惩罚项。

• 岭回归（Ridge Regression）：损失函数 =Σ(y_i - ŷ_i)^2 + λ * Σ(b_j^2)。惩罚项是系数的L2范数平方。它会让所有系数同时向零收缩，但不会将任何系数** exactly **压缩到零。适用于处理共线性。

  from sklearn.linear_model import Ridge model_ridge = Ridge(alpha=1.0) # alpha 就是 λ model_ridge.fit(X_train, y_train)

• Lasso回归（Lasso Regression）：损失函数 =Σ(y_i - ŷ_i)^2 + λ * Σ|b_j|。惩罚项是系数的L1范数。它倾向于产生稀疏解，即把一些不重要的变量的系数** exactly压缩到零，从而实现特征选择**。

  from sklearn.linear_model import Lasso model_lasso = Lasso(alpha=0.1) model_lasso.fit(X_train, y_train) # 查看哪些特征被筛掉了（系数为0） print(model_lasso.coef_)

参数λ（在sklearn中为alpha）控制惩罚力度，需要通过交叉验证来选取。

7.2 稳健回归（Robust Regression）

当数据中存在异常值（Outliers）时，OLS因为最小化平方和，会对异常值非常敏感（平方放大了大误差的影响）。稳健回归通过改变损失函数，降低异常值的权重。

• Huber回归：它对小误差使用平方损失，对大误差使用线性损失，从而减少异常值的影响。

  from sklearn.linear_model import HuberRegressor model_huber = HuberRegressor(epsilon=1.35) # epsilon是平方损失转向线性损失的阈值参数 model_huber.fit(X, y)

• RANSAC（随机抽样一致）：这是一种完全不同的思路。它反复随机抽取一个子集（假设是内点）来拟合模型，然后用这个模型去测试其他点，符合模型的点加入内点集。最终选择内点最多时拟合的模型。它对异常值有极强的鲁棒性。

  from sklearn.linear_model import RANSACRegressor from sklearn.linear_model import LinearRegression base_estimator = LinearRegression() model_ransac = RANSACRegressor(estimator=base_estimator, min_samples=0.5) model_ransac.fit(X, y) inlier_mask = model_ransac.inlier_mask_ # 标识哪些是内点 outlier_mask = ~inlier_mask

在实际项目中，我的经验是：永远先从OLS和残差分析开始。它简单、透明、易于解释。只有当诊断出明确的问题（如共线性、过拟合、异常值）时，再考虑引入这些更复杂的方法。记住，模型复杂度的提升，往往以牺牲可解释性为代价。