从数学原理到代码实现，理解回归的每一个细节
线性回归是最基本也最重要的统计学习方法。它假设因变量与自变量之间存在线性关系，用最小二乘法找到最佳拟合直线。
| 指标 | 公式 | 含义 | 理想值 |
|---|---|---|---|
| MSE | (1/n)Σ(ŷ-y)² | 均方误差 | 0 |
| RMSE | √MSE | 均方根误差（与y同量纲） | 0 |
| MAE | (1/n)Σ|ŷ-y| | 平均绝对误差 | 0 |
| R² | 1 - SS_res/SS_tot | 决定系数（解释方差比） | 1 |
| Adjusted R² | 1-(1-R²)(n-1)/(n-p-1) | 调整决定系数（惩罚特征数） | 1 |
梯度下降是机器学习最核心的优化算法。沿着损失函数梯度的反方向更新参数，逐步逼近最优解。
import numpy as np
np.random.seed(42)
# 生成数据: y = 2.5x + 3 + noise
n = 100
X = np.random.uniform(0, 10, n)
y = 2.5 * X + 3.0 + np.random.normal(0, 2, n)
# 标准化（梯度下降前必须做！）
X_mean, X_std = X.mean(), X.std()
y_mean, y_std = y.mean(), y.std()
X_norm = (X - X_mean) / X_std
y_norm = (y - y_mean) / y_std
# 梯度下降
w, b = 0.0, 0.0
lr = 0.01
losses = []
for epoch in range(1000):
y_pred = w * X_norm + b
error = y_pred - y_norm
dw = (2/n) * np.dot(error, X_norm)
db = (2/n) * np.sum(error)
w -= lr * dw
b -= lr * db
losses.append(np.mean(error**2))
# 反标准化
w_orig = w * (y_std / X_std)
b_orig = y_mean - w_orig * X_mean
print(f"学习参数: w={w_orig:.4f} (真实2.5), b={b_orig:.4f} (真实3.0)")
# R²计算
y_pred_orig = w_orig * X + b_orig
ss_res = np.sum((y - y_pred_orig)**2)
ss_tot = np.sum((y - y_mean)**2)
r_squared = 1 - ss_res / ss_tot
print(f"R²: {r_squared:.4f}") # 0.9403
print(f"MSE下降: {losses[0]:.4f} → {losses[-1]:.4f} (↓{(1-losses[-1]/losses[0])*100:.1f}%)")# 解析解 (Normal Equation): w = (XᵀX)⁻¹Xᵀy
X_with_bias = np.c_[np.ones(n), X] # 添加偏置列
w_analytic = np.linalg.inv(X_with_bias.T @ X_with_bias) @ X_with_bias.T @ y
print(f"解析解: w={w_analytic[1]:.4f}, b={w_analytic[0]:.4f}")
print(f"梯度下降: w={w_orig:.4f}, b={b_orig:.4f}")
print(f"差异: Δw={abs(w_analytic[1]-w_orig):.6f}")
# 性质对比
print("\n--- 解析解 vs 梯度下降 ---")
print("解析解: O(n³) 一次计算，小数据精确")
print("梯度下降: O(n·k·iter) 迭代，大数据可扩展")
print("当n>10万时，解析解矩阵求逆极慢，梯度下降更优")# 多元线性回归: y = w₁x₁ + w₂x₂ + b
np.random.seed(42)
n = 200
X_multi = np.random.randn(n, 3)
y_multi = 3 * X_multi[:, 0] - 2 * X_multi[:, 1] + 1.5 * X_multi[:, 2] + 5 + np.random.randn(n) * 0.5
# 解析解
X_bias = np.c_[np.ones(n), X_multi]
w_multi = np.linalg.inv(X_bias.T @ X_bias) @ X_bias.T @ y_multi
print(f"截距: {w_multi[0]:.4f}")
print(f"系数: {w_multi[1:]}")
# R²
y_pred_multi = X_bias @ w_multi
r2_multi = 1 - np.sum((y_multi - y_pred_multi)**2) / np.sum((y_multi - y_multi.mean())**2)
print(f"R²: {r2_multi:.4f}")
# Adjusted R²
p = X_multi.shape[1]
adj_r2 = 1 - (1 - r2_multi) * (n - 1) / (n - p - 1)
print(f"Adjusted R²: {adj_r2:.4f}")| 假设 | 检验方法 | 违反时的处理 |
|---|---|---|
| 线性关系 | 残差图（应随机分布） | 多项式回归/非线性变换 |
| 误差独立 | Durbin-Watson检验 | 时间序列模型 |
| 同方差性 | Breusch-Pagan检验 | 加权最小二乘/对数变换 |
| 正态性 | Shapiro-Wilk/Q-Q图 | 大样本下可放宽 |
| 无多重共线性 | VIF < 10 | 岭回归/Lasso/PCA |
# 多重共线性检测 (VIF)
from numpy.linalg import inv
def calc_vif(X, feature_idx):
"""计算第feature_idx个特征的VIF"""
y_vif = X[:, feature_idx]
X_others = np.delete(X, feature_idx, axis=1)
X_bias = np.c_[np.ones(len(y_vif)), X_others]
beta = inv(X_bias.T @ X_bias) @ X_bias.T @ y_vif
y_pred = X_bias @ beta
ss_res = np.sum((y_vif - y_pred)**2)
ss_tot = np.sum((y_vif - y_vif.mean())**2)
r2 = 1 - ss_res / ss_tot
return 1 / (1 - r2) if r2 < 1 else float('inf')
for i in range(X_multi.shape[1]):
vif = calc_vif(X_multi, i)
flag = '⚠️ 高共线性' if vif > 10 else '✅ OK'
print(f"特征{i} VIF: {vif:.2f} {flag}")| 方法 | 损失函数 | 特点 |
|---|---|---|
| OLS | MSE | 无正则化，可能过拟合 |
| Ridge (L2) | MSE + αΣwᵢ² | 收缩系数，不稀疏 |
| Lasso (L1) | MSE + αΣ|wᵢ| | 特征选择，稀疏解 |
| ElasticNet | MSE + α₁Σ|wᵢ| + α₂Σwᵢ² | 兼顾稀疏和分组效应 |
# Ridge回归手动实现
def ridge_regression(X, y, alpha=1.0):
"""w = (XᵀX + αI)⁻¹Xᵀy"""
X_bias = np.c_[np.ones(len(y)), X]
n_features = X_bias.shape[1]
I = np.eye(n_features)
I[0, 0] = 0 # 不正则化截距
w = inv(X_bias.T @ X_bias + alpha * I) @ X_bias.T @ y
return w
for alpha in [0, 0.1, 1, 10, 100]:
w_ridge = ridge_regression(X_multi, y_multi, alpha)
print(f"α={alpha:6.1f}: 截距={w_ridge[0]:.3f}, 系数={w_ridge[1:].round(3)}")