📈 回归基础

📖 回归基础：从数据到预测

回归分析是统计学和机器学习的基石。给定一组数据点，找到最优的直线（或超平面）来描述变量之间的关系，这就是回归。从手写最小二乘法到sklearn一行调用，理解回归原理是进入机器学习世界的第一步。

1. 手写最小二乘法(OLS)

import numpy as np

# 生成数据: y = 3x + 5 + 噪声
np.random.seed(42)
n = 100
X = np.random.uniform(0, 10, n)
y = 3 * X + 5 + np.random.normal(0, 2, n)

# 手写OLS: β = (X'X)^(-1) X'y
X_mat = np.column_stack([np.ones(n), X])  # 添加截距列
beta_hat = np.linalg.inv(X_mat.T @ X_mat) @ X_mat.T @ y
print(f"截距={beta_hat[0]:.4f}, 斜率={beta_hat[1]:.4f}")
# 截距≈5.0, 斜率≈3.0 → 完美恢复真实参数！

🤖 sklearn线性回归

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error

# 划分训练/测试集
X_train, X_test, y_train, y_test = train_test_split(
    X.reshape(-1,1), y, test_size=0.2, random_state=42)

# 训练模型
model = LinearRegression()
model.fit(X_train, y_train)

# 预测与评估
y_pred = model.predict(X_test)
print(f"截距={model.intercept_:.4f}, 斜率={model.coef_[0]:.4f}")
print(f"R²={r2_score(y_test, y_pred):.4f}")
print(f"RMSE={np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")
print(f"MAE={mean_absolute_error(y_test, y_pred):.4f}")

⚡ 梯度下降实现

当数据量大到矩阵求逆不可行时，梯度下降是替代方案：

def gradient_descent(X, y, lr=0.01, epochs=1000):
    m, n = X.shape
    w = np.zeros(n)
    b = 0
    for epoch in range(epochs):
        y_pred = X @ w + b
        dw = (2/m) * X.T @ (y_pred - y)
        db = (2/m) * np.sum(y_pred - y)
        w -= lr * dw
        b -= lr * db
        if (epoch+1) % 200 == 0:
            loss = np.mean((y_pred - y)**2)
            print(f"Epoch {epoch+1}: MSE={loss:.4f}")
    return w, b

# 多元回归: y = 2x₁ - 3x₂ + 0.5x₃ + 5
X_multi = np.random.randn(200, 3)
y_multi = 2*X_multi[:,0] - 3*X_multi[:,1] + 0.5*X_multi[:,2] + 5 + np.random.randn(200)*0.5

w_gd, b_gd = gradient_descent(X_train2, y_train2, lr=0.01, epochs=1000)
# 结果: w≈[1.967, -2.937, 0.559], b≈5.024 → 接近真实值！

📊 模型评估指标

残差分析

residuals = y_test - y_pred
print(f"残差均值: {residuals.mean():.6f}")  # 应≈0
print(f"残差标准差: {residuals.std():.4f}")

# 检查正态性
from scipy import stats
stat, p = stats.shapiro(residuals)

📐 2024-2025 回归前沿

• XGBoost/LightGBM：梯度提升回归，Kaggle常胜算法
• 神经网络回归：深度学习处理复杂非线性关系
• 分位数回归：预测条件分布的分位数，不只预测均值
• 贝叶斯回归：PyMC/Stan，参数的后验分布而非点估计
• 因果推断：Double ML，IV回归，分离因果效应

📊 回归诊断详解

import matplotlib.pyplot as plt
from scipy.stats import probplot

residuals = y_test - y_pred
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 残差 vs 预测值
axes[0].scatter(y_pred, residuals, alpha=0.6, color='#3b82f6')
axes[0].axhline(y=0, color='red', linestyle='--')

# 残差直方图
axes[1].hist(residuals, bins=20, color='#3b82f6', alpha=0.7)

# QQ图
probplot(residuals, plot=axes[2])

诊断检查清单

🧮 正则化回归

from sklearn.linear_model import Ridge, Lasso, ElasticNet
from sklearn.preprocessing import StandardScaler

# 标准化（正则化必须先标准化！）
scaler = StandardScaler()
X_s = scaler.fit_transform(X_train2)

# Ridge(L2): 系数衰减但不为零
ridge = Ridge(alpha=1.0)
ridge.fit(X_s, y_train2)

# Lasso(L1): 系数可变为零 → 自动特征选择
lasso = Lasso(alpha=0.1)
lasso.fit(X_s, y_train2)
print(f"非零特征: {(lasso.coef_!=0).sum()}")

# ElasticNet(L1+L2)
enet = ElasticNet(alpha=0.1, l1_ratio=0.5)
enet.fit(X_s, y_train2)

💻 完整实战代码

#!/usr/bin/env python3
# 回归基础 — 完整实战

import numpy as np
from numpy.linalg import lstsq
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error

# ============ 手写最小二乘法 ============
np.random.seed(42)
n = 100
X = np.random.uniform(0, 10, n)
y_true = 3 * X + 5
y = y_true + np.random.normal(0, 2, n)

X_mat = np.column_stack([np.ones(n), X])
beta_hat = np.linalg.inv(X_mat.T @ X_mat) @ X_mat.T @ y
print(f"手写OLS: 截距={beta_hat[0]:.4f}, 斜率={beta_hat[1]:.4f}")
print(f"理论值:  截距=5, 斜率=3")

y_pred_manual = X_mat @ beta_hat
ss_res = np.sum((y - y_pred_manual)**2)
ss_tot = np.sum((y - np.mean(y))**2)
r2_manual = 1 - ss_res / ss_tot
print(f"手写R²: {r2_manual:.4f}")

# ============ sklearn对比 ============
X_2d = X.reshape(-1, 1)
X_train, X_test, y_train, y_test = train_test_split(X_2d, y, test_size=0.2, random_state=42)
model = LinearRegression()
model.fit(X_train, y_train)
y_pred_sk = model.predict(X_test)
print(f"\nsklearn: 截距={model.intercept_:.4f}, 斜率={model.coef_[0]:.4f}")
print(f"sklearn R²: {r2_score(y_test, y_pred_sk):.4f}")
print(f"RMSE: {np.sqrt(mean_squared_error(y_test, y_pred_sk)):.4f}")
print(f"MAE: {mean_absolute_error(y_test, y_pred_sk):.4f}")

# ============ 多元回归 ============
np.random.seed(42)
X_multi = np.random.randn(200, 3)
y_multi = 2*X_multi[:,0] - 3*X_multi[:,1] + 0.5*X_multi[:,2] + 5 + np.random.randn(200)*0.5
X_train2, X_test2, y_train2, y_test2 = train_test_split(X_multi, y_multi, test_size=0.2, random_state=42)
model2 = LinearRegression()
model2.fit(X_train2, y_train2)
y_pred2 = model2.predict(X_test2)
print(f"\n多元回归系数: {model2.coef_.round(3)}")
print(f"截距: {model2.intercept_:.3f}")
print(f"多元R²: {r2_score(y_test2, y_pred2):.4f}")

# ============ 梯度下降 ============
def gradient_descent(X, y, lr=0.001, epochs=1000):
    m, n = X.shape
    w = np.zeros(n)
    b = 0
    for epoch in range(epochs):
        y_pred = X @ w + b
        dw = (2/m) * X.T @ (y_pred - y)
        db = (2/m) * np.sum(y_pred - y)
        w -= lr * dw
        b -= lr * db
        if (epoch+1) % 200 == 0:
            loss = np.mean((y_pred - y)**2)
            print(f"  Epoch {epoch+1}: MSE={loss:.4f}")
    return w, b

print(f"\n梯度下降(多元):")
w_gd, b_gd = gradient_descent(X_train2, y_train2, lr=0.01, epochs=1000)
y_pred_gd = X_test2 @ w_gd + b_gd
print(f"GD系数: {w_gd.round(3)}, 截距: {b_gd:.3f}")
print(f"GD R²: {r2_score(y_test2, y_pred_gd):.4f}")

print("\n✅ Python验证通过 — 最小二乘法手写+sklearn")

📝 课后练习

1. 实现Ridge和Lasso正则化回归
2. 用多项式特征处理非线性关系
3. 实现学习率衰减的梯度下降
4. 对真实数据做完整的回归分析（EDA+建模+诊断）
5. 实现贝叶斯线性回归（PyMC）

🗺️ 回归方法选择指南

正则化方法对比