AI 工程基础体系 · 第 42/100 篇。内容覆盖机器学习、深度学习与生成式 AI;模型、数据、评测、权限和成本会作为同一生产系统处理。

线性回归完整原理:最小二乘、正则化、诊断与置信区间

线性回归(linear regression)用一个关于输入特征的线性函数预测连续数值。给定样本 ii 的特征向量 xiRpx_i\in\mathbb{R}^p 和目标值 yiRy_i\in\mathbb{R},最基本的模型是:

yi=β0+β1xi1++βpxip+εiy_i=\beta_0+\beta_1x_{i1}+\cdots+\beta_px_{ip}+\varepsilon_i

其中:

  • β0\beta_0 是截距;
  • βj\beta_j 是第 jj 个特征的系数;
  • εi\varepsilon_i 是模型未解释的误差;
  • “线性”指对参数 β\beta 线性,而不一定要求原始变量只能以一次方形式出现。

例如:

y=β0+β1x+β2x2+εy=\beta_0+\beta_1x+\beta_2x^2+\varepsilon

对参数 β0,β1,β2\beta_0,\beta_1,\beta_2 仍然是线性回归,只是把 x2x^2 作为了一个新的特征。相反:

y=β0eβ1x+εy=\beta_0e^{\beta_1x}+\varepsilon

通常称为非线性参数模型,因为参数出现在指数函数中。


一、从数据到矩阵形式

nn 个样本写成矩阵:

X=[1x11x1p1x21x2p1xn1xnp],y=[y1y2yn]X= \begin{bmatrix} 1 & x_{11} & \cdots & x_{1p}\\ 1 & x_{21} & \cdots & x_{2p}\\ \vdots & \vdots & & \vdots\\ 1 & x_{n1} & \cdots & x_{np} \end{bmatrix}, \quad y= \begin{bmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{bmatrix}

第一列全为 11,用于表示截距。令:

θ=[β0β1βp]\theta= \begin{bmatrix} \beta_0\\ \beta_1\\ \vdots\\ \beta_p \end{bmatrix}

则模型可以简写为:

y=Xθ+εy=X\theta+\varepsilon

对单个样本,预测值是:

y^i=xiθ\hat y_i=x_i^\top\theta

这里的 xix_i 已经包含了第一维常数 11

线性回归解决的具体问题

训练阶段要从数据中选择一个参数向量 θ\theta,使预测值 y^=Xθ\hat y=X\theta 尽可能接近真实值 yy。最常用的距离是平方误差:

SSE(θ)=i=1n(yiy^i)2=yXθ22\operatorname{SSE}(\theta) =\sum_{i=1}^n(y_i-\hat y_i)^2 =\|y-X\theta\|_2^2

其中:

  • SSE\operatorname{SSE} 是残差平方和;
  • r=yXθr=y-X\theta 是残差向量;
  • 2\|\cdot\|_2 是欧氏范数。

平方误差会放大较大的错误。例如误差为 11 时平方为 11,误差为 55 时平方为 2525。因此,最小二乘会特别重视大残差,这既是它有效的原因,也是它容易受异常值影响的原因。


二、最小二乘估计的推导

2.1 正规方程

定义损失函数:

L(θ)=yXθ22=(yXθ)(yXθ)L(\theta)=\|y-X\theta\|_2^2 =(y-X\theta)^\top(y-X\theta)

展开:

L(θ)=yy2θXy+θXXθL(\theta) =y^\top y-2\theta^\top X^\top y+\theta^\top X^\top X\theta

θ\theta 求梯度:

θL=2Xy+2XXθ\nabla_\theta L =-2X^\top y+2X^\top X\theta

令梯度为零:

XXθ^=XyX^\top X\hat\theta=X^\top y

这组方程称为正规方程

如果 XX 的列满秩,则 XXX^\top X 可逆,得到:

θ^=(XX)1Xy\boxed{\hat\theta=(X^\top X)^{-1}X^\top y}

这就是普通最小二乘(ordinary least squares,OLS)估计量。

2.2 为什么是最小值而不是最大值

平方损失的 Hessian 矩阵为:

2L=2XX\nabla^2L=2X^\top X

对任意向量 vv

vXXv=Xv220v^\top X^\top Xv=\|Xv\|_2^2\ge 0

因此 XXX^\top X 是半正定矩阵,损失函数是凸函数。任何满足正规方程的点都是全局最小值。

如果 XX 列不满秩,参数解可能不唯一,但预测值 Xθ^X\hat\theta 仍可能唯一。实际实现通常使用 QR 分解或 SVD,而不是直接计算 (XX)1(X^\top X)^{-1},因为显式求逆数值稳定性较差。


三、几何解释:投影为何产生正交残差

XθX\theta 的所有可能取值组成 XX 的列空间:

C(X)={Xθ:θRp+1}\mathcal{C}(X)=\{X\theta:\theta\in\mathbb{R}^{p+1}\}

最小二乘寻找列空间中距离 yy 最近的点:

y^=Xθ^\hat y=X\hat\theta

因此,y^\hat yyyC(X)\mathcal{C}(X) 的正交投影。投影误差 r=yy^r=y-\hat y 必须与列空间正交:

Xr=0X^\top r=0

这正好就是正规方程:

X(yXθ^)=0X^\top(y-X\hat\theta)=0

如果包含截距,第一列是全 1,因此:

1r=0\mathbf{1}^\top r=0

也就是说,包含截距的 OLS 拟合中,残差和为零,训练集上的残差均值为零。这个性质只适用于相应条件成立的训练拟合,不能推出测试集残差均值也为零。


四、一个完整的最小二乘算例

考虑三个观测点:

(x,y)=(1,2),(2,3),(3,5)(x,y)=(1,2),(2,3),(3,5)

模型为:

y=β0+β1x+εy=\beta_0+\beta_1x+\varepsilon

写成矩阵:

X=[111213],y=[235]X= \begin{bmatrix} 1&1\\ 1&2\\ 1&3 \end{bmatrix}, \quad y= \begin{bmatrix} 2\\3\\5 \end{bmatrix}

先计算:

XX=[36614],Xy=[1023]X^\top X= \begin{bmatrix} 3&6\\ 6&14 \end{bmatrix}, \quad X^\top y= \begin{bmatrix} 10\\ 23 \end{bmatrix}

正规方程为:

[36614][β^0β^1]=[1023]\begin{bmatrix} 3&6\\ 6&14 \end{bmatrix} \begin{bmatrix} \hat\beta_0\\ \hat\beta_1 \end{bmatrix} = \begin{bmatrix} 10\\ 23 \end{bmatrix}

解得:

β^1=1.5,β^0=13\hat\beta_1=1.5,\qquad \hat\beta_0=\frac13

因此:

y^=13+1.5x\hat y=\frac13+1.5x

三个预测值为:

y^=[1.83333.33334.8333]\hat y= \begin{bmatrix} 1.8333\\ 3.3333\\ 4.8333 \end{bmatrix}

残差为:

r=yy^=[0.16670.33330.1667]r=y-\hat y= \begin{bmatrix} 0.1667\\ -0.3333\\ 0.1667 \end{bmatrix}

残差平方和:

SSE=0.16672+(0.3333)2+0.166720.1667\operatorname{SSE} =0.1667^2+(-0.3333)^2+0.1667^2 \approx0.1667

残差和为零,且:

ixiri0\sum_i x_ir_i \approx 0

说明残差同时与截距列和特征列正交。

这个算例也暴露了自由度问题:估计了两个参数,却只有三个观测,因此残差自由度为:

ν=n(p+1)=32=1\nu=n-(p+1)=3-2=1

样本太少时,系数和置信区间都非常不稳定。能够计算不代表统计结论可靠。


五、系数的含义与“控制其他变量不变”

在模型:

y^=β^0+β^1x1++β^pxp\hat y=\hat\beta_0+\hat\beta_1x_1+\cdots+\hat\beta_px_p

中,β^j\hat\beta_j 表示:在其他输入特征保持不变时,xjx_j 增加一个单位,预测值平均改变 β^j\hat\beta_j 个单位。

“保持其他变量不变”是数学模型中的条件解释,不等同于现实世界中的因果结论。要声称“提高 xjx_j 会导致 yy 增加”,还需要处理混杂变量、反向因果、选择偏差和实验设计等问题。

如果对特征做了标准化:

zj=xjxˉjsjz_j=\frac{x_j-\bar x_j}{s_j}

则对应系数表示 xjx_j 增加一个标准差时预测值的变化。标准化有利于比较系数大小,也有利于正则化,但会改变系数的单位,不能直接与原始尺度上的系数混用。


六、平方损失与概率模型的关系

假设误差独立同分布,并服从正态分布:

εiN(0,σ2)\varepsilon_i\sim\mathcal N(0,\sigma^2)

则:

yixiN(xiθ,σ2)y_i\mid x_i\sim\mathcal N(x_i^\top\theta,\sigma^2)

对数似然为:

logp(yX,θ,σ2)=n2log(2πσ2)12σ2yXθ22\log p(y\mid X,\theta,\sigma^2) =-\frac{n}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2}\|y-X\theta\|_2^2

σ2\sigma^2 固定时,最大化似然等价于最小化 SSE。因此,OLS 既可以理解为几何上的投影,也可以理解为高斯噪声假设下的最大似然估计。

但要注意:

  • 高斯误差不是 OLS 能计算的必要条件;
  • 正态性主要用于精确的小样本检验和置信区间;
  • 在大样本下,某些条件下可以使用渐近近似;
  • 如果目标值是类别、计数或概率,通常应考虑逻辑回归、泊松回归、Beta 回归等广义线性模型,而不是强行使用普通线性回归。

七、OLS 的统计假设及其后果

7.1 线性条件均值

关键条件通常写成:

E[εX]=0\mathbb E[\varepsilon\mid X]=0

它表示在给定特征后,误差的平均值为零。该条件保证:

E[θ^X]=θ\mathbb E[\hat\theta\mid X]=\theta

也就是 OLS 对参数无偏。

如果遗漏了与输入相关的重要变量,或者特征与误差相关,则可能有:

E[εX]0\mathbb E[\varepsilon\mid X]\ne0

此时增加样本并不能自动消除偏差,这称为遗漏变量偏差或内生性问题。

7.2 独立性

如果不同样本的误差相关,例如时间序列中的相邻观测,普通标准误可能错误。系数点估计仍可能有意义,但置信区间和显著性检验会失真。

随机划分时间序列数据还可能造成未来信息泄露。时间数据通常应按时间切分,并考虑滞后特征、滚动验证或时间序列误差结构。

7.3 同方差性

同方差要求:

Var(εiX)=σ2\operatorname{Var}(\varepsilon_i\mid X)=\sigma^2

如果误差方差随输入改变,就是异方差。例如收入越高,消费金额的绝对波动可能越大。

异方差通常不必然破坏 OLS 系数的无偏性,但会使经典标准误、t 检验和置信区间不可靠。可采用:

  • 异方差稳健标准误;
  • 加权最小二乘;
  • 对目标做对数或其他方差稳定变换;
  • 显式建模条件方差。

加权最小二乘的目标为:

minθiwi(yixiθ)2\min_\theta\sum_iw_i(y_i-x_i^\top\theta)^2

矩阵解为:

θ^WLS=(XWX)1XWy\hat\theta_{\mathrm{WLS}} =(X^\top WX)^{-1}X^\top Wy

其中 WW 通常是对角矩阵。若 wiw_i 与样本误差方差的倒数成比例,估计会更有效;权重设置错误则可能引入新的问题。

7.4 多重共线性

如果两个特征高度相关,XXX^\top X 会接近奇异。此时:

  • 系数方差增大;
  • 系数符号可能反复变化;
  • 单个系数的解释变得不稳定;
  • 整体预测可能仍然准确。

多重共线性不是简单的“模型错误”。如果目标是预测,可以接受一定共线性;如果目标是解释单个变量,则需要合并特征、删除冗余特征、使用主成分或正则化,并明确解释对象已经改变。


八、R2R^2、MSE 与评测边界

常见训练指标包括:

MSE=1ni(yiy^i)2\operatorname{MSE} =\frac1n\sum_i(y_i-\hat y_i)^2

RMSE=MSE\operatorname{RMSE}=\sqrt{\operatorname{MSE}}

RMSE 与目标值同单位,更容易解释。

带截距模型常用决定系数:

R2=1i(yiy^i)2i(yiyˉ)2R^2 =1-\frac{\sum_i(y_i-\hat y_i)^2} {\sum_i(y_i-\bar y)^2}

其中分母是总平方和。R2R^2 表示模型相对“始终预测训练集均值”的基线减少了多少平方误差。

边界包括:

  1. 测试集上的 R2R^2 可以为负,表示还不如均值基线。
  2. 增加特征不会降低训练集 R2R^2,即使新特征只是噪声。
  3. R2R^2 不证明因果关系,也不证明未来数据表现好。
  4. 不同目标尺度或不同数据集之间,RMSE 不能直接比较。
  5. 对异常值敏感的平方损失,可能掩盖大多数样本上的典型表现。

生产评测应在独立验证数据上进行,并固定数据切分、目标定义、时间窗口和缺失值处理方式。对于高成本预测错误,还应同时报告分位数误差、业务加权损失或分段误差。


九、正则化:为什么需要它

当特征很多、样本较少、特征高度相关或存在过拟合风险时,单纯最小化训练误差可能得到方差很大的系数。

正则化在训练目标中加入对参数大小的惩罚:

训练损失+λ复杂度惩罚\text{训练损失}+\lambda\cdot\text{复杂度惩罚}

λ0\lambda\ge0 控制惩罚强度:

  • λ=0\lambda=0 时退化为 OLS;
  • λ\lambda 增大时,参数通常更小;
  • λ\lambda 太大则可能产生明显欠拟合。

正则化通常改变的是偏差—方差权衡:允许少量偏差,换取更低的估计方差。

9.1 Ridge 回归

Ridge 使用 L2 惩罚:

θ^ridge=argminθ{yXθ22+λβ22}\hat\theta_{\text{ridge}} = \arg\min_\theta \left\{ \|y-X\theta\|_2^2+\lambda\|\beta\|_2^2 \right\}

其中通常不惩罚截距 β0\beta_0β\beta 表示不含截距的系数。其闭式解为:

β^ridge=(XX+λI)1Xy\hat\beta_{\text{ridge}} =(X^\top X+\lambda I)^{-1}X^\top y

如果显式区分截距,惩罚矩阵应为:

P=[000I]P= \begin{bmatrix} 0&0\\ 0&I \end{bmatrix}

解为:

θ^=(XX+λP)1Xy\hat\theta=(X^\top X+\lambda P)^{-1}X^\top y

Ridge 的重要作用是让矩阵更稳定。即使 XXX^\top X 接近奇异,增加 λI\lambda I 后也通常更容易求解。Ridge 会把相关特征的系数一起压小,但一般不会精确变成零。

9.2 Lasso 回归

Lasso 使用 L1 惩罚:

θ^lasso=argminθ{yXθ22+λβ1}\hat\theta_{\text{lasso}} = \arg\min_\theta \left\{ \|y-X\theta\|_2^2+\lambda\|\beta\|_1 \right\}

其中:

β1=jβj\|\beta\|_1=\sum_j|\beta_j|

L1 惩罚在零点不可导,产生“尖角”,因此一些系数会精确为零,具有稀疏选择效果。

但 Lasso 的选择并不等价于发现真正的因果变量。高度相关的一组特征中,它可能任意保留其中一个,或者随着数据变化在它们之间切换。

9.3 Elastic Net

Elastic Net 同时使用 L1 和 L2:

minθ{yXθ22+λ[αβ1+(1α)β22]}\min_\theta \left\{ \|y-X\theta\|_2^2+ \lambda\left[ \alpha\|\beta\|_1+ (1-\alpha)\|\beta\|_2^2 \right] \right\}

其中 α[0,1]\alpha\in[0,1]

  • α=1\alpha=1 接近 Lasso;
  • α=0\alpha=0 接近 Ridge;
  • 中间值兼顾稀疏性和对相关特征的稳定处理。

不同库对损失是否除以 nn、惩罚系数如何缩放的约定可能不同,因此不能直接把一个实现中的 λ\lambda 数值照搬到另一个实现。

9.4 为什么必须缩放特征

假设一个特征单位是元,另一个单位是平方米。L1/L2 惩罚作用于系数,而不是特征对预测的实际贡献。若特征尺度差异很大,系数大小不能公平比较,正则化会偏向惩罚某些特征。

常见流程是:

  1. 用训练集估计均值和标准差;
  2. 对训练、验证、测试数据使用同一组统计量;
  3. 在缩放后的特征上训练正则化模型。

不能先用全量数据计算均值和标准差,否则测试集信息会泄露到训练流程。


十、正则化参数的选择

λ\lambda 不能根据测试集表现反复挑选。常见做法是训练集内部交叉验证:

  1. 先划分训练集和最终测试集;
  2. 在训练集内部进行 K 折交叉验证;
  3. 对候选 λ\lambda 训练模型;
  4. 选择验证损失最小或符合业务约束的值;
  5. 用选定的超参数在整个训练集上重新拟合;
  6. 最后只在测试集上评估一次。

如果数据有时间顺序、用户分组或实体重复,普通随机 K 折可能把同一实体或未来信息同时放入训练和验证,导致过于乐观的结果。应采用时间切分、GroupKFold 或其他符合数据生成过程的切分方式。


十一、从系数到置信区间

置信区间描述的是一个统计程序在重复抽样下覆盖未知参数的频率性质,不是说“这个固定参数以某个概率位于当前区间内”。

11.1 OLS 系数的分布

在固定设计矩阵 XX、模型正确且:

E[εX]=0,Var(εX)=σ2I\mathbb E[\varepsilon\mid X]=0,\qquad \operatorname{Var}(\varepsilon\mid X)=\sigma^2I

时:

θ^=(XX)1Xy\hat\theta=(X^\top X)^{-1}X^\top y

代入 y=Xθ+εy=X\theta+\varepsilon

θ^=θ+(XX)1Xε\hat\theta =\theta+(X^\top X)^{-1}X^\top\varepsilon

因此:

E[θ^X]=θ\mathbb E[\hat\theta\mid X]=\theta

并且:

Var(θ^X)=σ2(XX)1\operatorname{Var}(\hat\theta\mid X) =\sigma^2(X^\top X)^{-1}

如果误差进一步服从正态分布,则:

θ^XN(θ,σ2(XX)1)\hat\theta\mid X \sim \mathcal N\left( \theta,\, \sigma^2(X^\top X)^{-1} \right)

11.2 估计噪声方差

真实的 σ2\sigma^2 通常未知,用残差估计:

σ^2=SSEnk\hat\sigma^2 =\frac{\operatorname{SSE}}{n-k}

其中 k=p+1k=p+1 是包含截距的参数数量。因此第 jj 个系数的标准误为:

SE(θ^j)=σ^2[(XX)1]jj\operatorname{SE}(\hat\theta_j) = \sqrt{ \hat\sigma^2 \left[(X^\top X)^{-1}\right]_{jj} }

分母使用 nkn-k,而不是 nn,因为参数估计消耗了自由度。

11.3 系数置信区间

在正态误差假设下,系数的 100(1α)%100(1-\alpha)\% 置信区间为:

θ^j±t1α/2,nkSE(θ^j)\hat\theta_j \pm t_{1-\alpha/2,n-k} \operatorname{SE}(\hat\theta_j)

其中 t1α/2,nkt_{1-\alpha/2,n-k} 是自由度为 nkn-k 的 t 分布分位数。

区间包含零时,不能简单说“该变量没有作用”。更准确的表述是:在当前模型、样本和误差假设下,数据没有提供足够证据拒绝系数为零。区间宽度还受到样本量、噪声、共线性和特征尺度影响。


十二、均值响应区间与新样本预测区间

这两个区间经常被混淆。

对一个新的特征向量 x0x_0,其平均响应为:

μ0=x0θ\mu_0=x_0^\top\theta

估计的平均响应为:

μ^0=x0θ^\hat\mu_0=x_0^\top\hat\theta

其方差为:

Var(μ^0X)=σ2x0(XX)1x0\operatorname{Var}(\hat\mu_0\mid X) = \sigma^2x_0^\top(X^\top X)^{-1}x_0

因此均值响应的区间为:

μ^0±t1α/2,nkσ^x0(XX)1x0\hat\mu_0 \pm t_{1-\alpha/2,n-k} \hat\sigma \sqrt{x_0^\top(X^\top X)^{-1}x_0}

但新观测还包含一个新的随机误差:

y0=x0θ+ε0y_0=x_0^\top\theta+\varepsilon_0

所以新样本预测误差的方差是:

Var(y0μ^0X)=σ2[1+x0(XX)1x0]\operatorname{Var}(y_0-\hat\mu_0\mid X) = \sigma^2 \left[ 1+x_0^\top(X^\top X)^{-1}x_0 \right]

预测区间为:

y^0±t1α/2,nkσ^1+x0(XX)1x0\boxed{ \hat y_0 \pm t_{1-\alpha/2,n-k} \hat\sigma \sqrt{ 1+x_0^\top(X^\top X)^{-1}x_0 } }

其中额外的 11 就是新样本自身噪声造成的。因此预测区间一定比同置信水平的均值响应区间更宽。样本数量增加时,参数不确定性会下降,但不可约的观测噪声不会自动消失。


十三、线性回归诊断:先看残差,而不是只看分数

诊断的目标是发现模型假设与数据之间的冲突。一个高 R2R^2 模型仍可能存在严重异方差、泄露、异常点或外推风险。

13.1 残差—拟合值图

绘制:

ri=yiy^ir_i=y_i-\hat y_i

相对于 y^i\hat y_i 的散点图。

理想情况下,点应围绕零随机散布。常见模式及含义:

  • U 形或倒 U 形:可能遗漏非线性;
  • 漏斗形:可能存在异方差;
  • 连续上升或下降:可能有时间依赖或系统性偏差;
  • 分组条带:可能遗漏类别变量或分层结构。

若关系明显非线性,可以增加多项式、样条、交互项,或改用树模型、广义加性模型等。增加高次多项式虽然能降低训练误差,却可能在边界剧烈振荡,因此必须依靠验证集和正则化控制复杂度。

13.2 杠杆值

帽子矩阵为:

H=X(XX)1XH=X(X^\top X)^{-1}X^\top

拟合值满足:

y^=Hy\hat y=Hy

ii 个对角元素 hiih_{ii} 称为杠杆值,反映样本特征向量在 XX 空间中是否远离其他样本。高杠杆点不一定是错误数据,但它有能力显著影响拟合直线。

在包含 kk 个参数时,杠杆值平均为:

1nihii=kn\frac{1}{n}\sum_i h_{ii}=\frac{k}{n}

经验上会关注明显高于平均水平的点,但不存在只凭一个固定阈值就能判定异常的普适规则。

13.3 学生化残差与 Cook 距离

残差大并不一定意味着影响大;高杠杆点即使残差不大,也可能改变模型。Cook 距离同时考虑残差和杠杆值,用于衡量删除某个样本后拟合结果的变化。

诊断流程通常是:

  1. 找出大残差样本;
  2. 找出高杠杆样本;
  3. 找出 Cook 距离较大的样本;
  4. 回查原始数据、单位、采集时间和标签;
  5. 分别比较保留、修正、剔除或使用稳健方法后的结果。

不能因为一个点让指标变差就直接删除。它可能是录入错误,也可能是真实而重要的极端业务场景。

13.4 正态性检查的正确用途

Q-Q 图可以检查残差是否近似正态。残差正态性对 OLS 点估计不是必要条件,但对小样本下的精确 t 区间和检验有帮助。

“残差看起来不正态”不自动意味着模型不能用于预测。若预测性能稳定,可以继续使用;但小样本置信区间、尾部风险和异常值敏感性需要更谨慎,必要时可使用变换、稳健标准误、Bootstrap 或稳健回归。


十四、异常值、稳健回归与替代损失

OLS 最小化平方误差,对大残差非常敏感。若误差分布有重尾或数据中存在少量异常值,可考虑:

Huber 损失

Huber 损失在小误差时近似平方损失,在大误差时近似绝对值损失:

Lδ(r)={12r2,rδδ(r12δ),r>δL_\delta(r)= \begin{cases} \frac12r^2,&|r|\le\delta\\ \delta(|r|-\frac12\delta),&|r|>\delta \end{cases}

它减少极端残差对参数的影响,但不等价于自动识别并删除异常样本。

Quantile 回归

分位数回归估计条件分位数,而不是条件均值。其 pinball 损失为:

ρτ(r)={τr,r0(τ1)r,r<0\rho_\tau(r)= \begin{cases} \tau r,&r\ge0\\ (\tau-1)r,&r<0 \end{cases}

例如 τ=0.5\tau=0.5 对应条件中位数,适用于关注上界、下界或不对称风险的场景。

改变损失函数后,系数含义、评测指标和区间解释也会改变,不能仍然套用 OLS 的标准误公式。


十五、可执行示例:拟合、系数区间和预测区间

下面示例使用 NumPy、scikit-learn 和 SciPy,演示一个含截距的一元线性回归,并手动计算系数置信区间、均值响应区间和预测区间。

安装依赖:

python -m pip install numpy scipy scikit-learn

代码:

import numpy as np
from scipy.stats import t
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

# 三个训练样本
x = np.array([1.0, 2.0, 3.0])
y = np.array([2.0, 3.0, 5.0])
X = x.reshape(-1, 1)

# fit_intercept=True 表示估计截距
model = LinearRegression(fit_intercept=True)
model.fit(X, y)

y_hat = model.predict(X)
residuals = y - y_hat

n = len(y)
k = X.shape[1] + 1  # 1 个特征 + 1 个截距
dof = n - k

# 加入截距列,构造设计矩阵
X_design = np.column_stack([np.ones(n), x])

sse = np.sum(residuals ** 2)
sigma2_hat = sse / dof
cov_theta = sigma2_hat * np.linalg.inv(X_design.T @ X_design)
se_theta = np.sqrt(np.diag(cov_theta))

alpha = 0.05
t_critical = t.ppf(1 - alpha / 2, dof)

theta_hat = np.array([model.intercept_, model.coef_[0]])
ci_theta = np.column_stack([
    theta_hat - t_critical * se_theta,
    theta_hat + t_critical * se_theta
])

# 对 x0=2.5 计算区间
x0 = np.array([1.0, 2.5])  # 第一维是截距
y0_hat = x0 @ theta_hat
leverage_term = x0 @ np.linalg.inv(X_design.T @ X_design) @ x0

mean_se = np.sqrt(sigma2_hat * leverage_term)
prediction_se = np.sqrt(sigma2_hat * (1 + leverage_term))

mean_interval = np.array([
    y0_hat - t_critical * mean_se,
    y0_hat + t_critical * mean_se
])

prediction_interval = np.array([
    y0_hat - t_critical * prediction_se,
    y0_hat + t_critical * prediction_se
])

print("intercept =", model.intercept_)
print("coef      =", model.coef_[0])
print("predicted =", y_hat)
print("SSE       =", sse)
print("RMSE      =", mean_squared_error(y, y_hat, squared=False))
print("R2        =", r2_score(y, y_hat))
print("coefficient CI [intercept, coef] =")
print(ci_theta)
print("x0=2.5 fitted mean =", y0_hat)
print("mean response CI   =", mean_interval)
print("new observation PI =", prediction_interval)

在当前数据上,关键结果应接近:

intercept = 0.3333333333
coef      = 1.5
predicted = [1.83333333 3.33333333 4.83333333]
SSE       = 0.1666666667

这个示例的置信区间会非常宽,因为自由度只有 11。这不是代码错误,而是样本信息不足的直接结果。若把训练集上的区间误当成稳定的生产不确定性,就会得到过度自信的结论。

代码中使用 np.linalg.inv 是为了展示公式。生产实现不应在大矩阵上直接显式求逆;应优先使用库提供的数值稳定求解器、QR 分解或 SVD。对于更复杂的诊断和稳健标准误,应使用经过验证的统计工具,而不是手动复制部分公式。


十六、scikit-learn 中的生产训练管线

对于包含缩放和正则化的模型,预处理必须和模型绑定,避免交叉验证或部署时发生不一致:

import numpy as np

from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split, GridSearchCV, KFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error, r2_score

data = load_diabetes()
X, y = data.data, data.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

pipeline = Pipeline([
    ("scale", StandardScaler()),
    ("model", Ridge())
])

# 参数搜索只在训练集内部完成
search = GridSearchCV(
    estimator=pipeline,
    param_grid={"model__alpha": np.logspace(-4, 4, 30)},
    scoring="neg_root_mean_squared_error",
    cv=KFold(n_splits=5, shuffle=True, random_state=42),
    n_jobs=-1
)

search.fit(X_train, y_train)

pred = search.predict(X_test)

print("best alpha =", search.best_params_["model__alpha"])
print("test RMSE  =", mean_squared_error(y_test, pred, squared=False))
print("test R2    =", r2_score(y_test, pred))

这里的因果关系是:

  • Pipeline 让每个交叉验证折只用该折的训练部分拟合 StandardScaler
  • GridSearchCV 只用训练集选择 alpha
  • 测试集直到最后才参与评估;
  • Ridgealpha 是实现中的正则化强度参数,具体目标函数缩放应以当前 scikit-learn 版本文档为准,不能与其他库的 λ\lambda 数值直接比较。

如果在交叉验证前先对全部数据调用 fit_transform,测试折的均值和标准差就会泄露到训练折,导致验证结果偏乐观。


十七、正则化模型的置信区间为何不能直接套用 OLS 公式

OLS 的经典区间依赖于:

θ^=(XX)1Xy\hat\theta=(X^\top X)^{-1}X^\top y

而 Ridge 和 Lasso 的估计量不同。

Ridge 的估计有偏:

E[β^ridgeX]=(XX+λI)1XXβ\mathbb E[\hat\beta_{\text{ridge}}\mid X] = (X^\top X+\lambda I)^{-1}X^\top X\beta

它通常降低方差,但不再满足 OLS 的无偏公式。Lasso 还涉及变量选择,选择事件本身会影响推断分布。

因此,对正则化模型不能简单地:

  1. 取正则化后的系数;
  2. 使用 σ^2(XX)1\hat\sigma^2(X^\top X)^{-1}
  3. 宣称得到标准 OLS 置信区间。

如果确实需要不确定性估计,应明确目标和方法,例如:

  • Bootstrap;
  • 贝叶斯线性模型;
  • 专门的选择后推断方法;
  • 通过重复重采样评估系数和预测分布;
  • 只报告经过验证的预测区间,而不把系数当作无偏解释量。

Bootstrap 也不是无条件可靠:时间序列、分组数据、强依赖数据和极端小样本需要采用相应的重采样方案。


十八、外推、分布偏移与数据泄露

线性回归在训练数据覆盖范围内进行插值时,通常比在范围外外推更可信。若训练数据中的年龄范围是 20 到 60 岁,却用模型预测 150 岁对象,直线会继续延伸,但数据并没有支持这种行为。

生产系统还要检查:

  • 训练和线上特征的单位、编码和缺失规则是否一致;
  • 特征是否在预测时刻真实可用;
  • 是否误用了未来字段、人工审核结果或目标发生后的统计量;
  • 新用户、新地区、新设备是否超出训练分布;
  • 目标定义是否随时间变化;
  • 训练数据、模型文件和预测日志是否含有受限个人信息。

数据权限不是模型拟合之后才处理的问题。谁可以读取训练数据、谁可以查看带真实标签的评测结果、谁可以下载模型和预测日志,都会影响系统风险。线性模型参数少、推理成本低,但低成本不代表可以跳过数据治理和访问控制。


十九、线性回归在深度学习与生成式 AI 系统中的位置

线性回归不是深度学习的对立面。神经网络最后一层常见的回归头就是:

y^=Wh+b\hat y=W h+b

其中 hh 是前面网络产生的表示,W,bW,b 是线性层参数。如果使用均方误差训练,这一部分仍遵循平方损失的基本机制,只是输入表示 hh 由深度网络学习得到。

在线性回归中,特征工程决定了模型能表达的函数族;在深度模型中,隐藏层承担了部分表示学习功能。但以下问题依然存在:

  • 训练目标与线上目标是否一致;
  • 误差是否受异常值影响;
  • 评测集是否泄露;
  • 置信区间或预测区间是否校准;
  • 模型输出是否被错误解释为因果结论;
  • 推理成本、延迟、权限和日志保留是否满足生产要求。

在生成式 AI 系统中,线性回归还可能作为成本、延迟、质量分数或资源使用量的预测器。例如用请求特征预测推理成本时,目标值往往异方差、重尾且受模型版本影响,此时应按版本和时间切分数据,并检查预测区间在不同流量段是否仍然覆盖,而不能只看总体 RMSE。


二十、常见误解与对应失败表现

误解一:线性回归必须画成直线

错误。线性指参数线性。加入 x2x^2、交互项、对数变换后的特征后,仍可能使用线性回归估计参数。

误解二:系数显著就代表因果关系

错误。显著性只能说明在特定模型和假设下,数据与某个参数值存在统计差异,不能自动排除混杂和反向因果。

误解三:正则化只是让训练更慢

错误。正则化改变了优化目标和最终估计量,可能增加训练误差,却降低测试误差和系数方差。

误解四:删除所有大残差点

错误。大残差可能是真实稀有事件。直接删除会让模型只适用于“正常样本”,并且可能造成选择偏差。

误解五:训练集 R2R^2 很高就说明模型可上线

错误。模型可能记住了噪声,或者使用了线上不可用的未来信息。必须在符合实际生成过程的独立数据上验证。

误解六:95% 置信区间表示新样本有 95% 概率落入其中

错误。系数置信区间描述参数推断;新样本应使用预测区间。即使是预测区间,也只有在模型和覆盖假设合理时才具有相应频率解释。


二十一、一个可靠的分析顺序

对于一个新的连续值预测问题,可以按以下因果顺序进行:

  1. 定义目标:明确预测时点、单位、时间范围和允许的误差。
  2. 检查数据生成过程:确认样本是否独立,是否存在时间、用户或组结构。
  3. 建立基线:至少比较均值预测和简单业务规则。
  4. 拟合 OLS:先理解特征、残差和系数,不要一开始就隐藏问题。
  5. 绘制诊断图:检查非线性、异方差、异常值、杠杆点和时间依赖。
  6. 处理结构性问题:增加合理变换、交互项、权重或稳健损失。
  7. 使用正则化:通过训练集内部交叉验证选择 Ridge、Lasso 或 Elastic Net 强度。
  8. 计算不确定性:区分系数区间、均值响应区间和新观测预测区间。
  9. 进行独立评测:使用不泄露的时间、实体或分组切分。
  10. 上线后监控:同时监控输入分布、残差、覆盖率、成本、延迟和权限事件。

线性回归的价值不只在于得到一条拟合直线。它把预测、优化、误差结构、参数不确定性和模型诊断放在同一个可计算框架中。掌握最小二乘的几何与概率解释,理解正则化如何改变偏差—方差权衡,并能区分诊断信号与业务因果,才是将线性模型可靠地用于机器学习和生产系统的基础。


系列导航与关联阅读

官方资料

本文依据研究论文、标准组织与主流框架官方文档重新梳理;正文、示例与工程清单由 WR BLOG 编写。