从零实现线性回归：波士顿房价预测实战指南

在机器学习领域，线性回归堪称算法界的"Hello World"。这个看似简单的模型，却蕴含着丰富的数学原理和实践技巧。本文将带你用Python从零开始构建线性回归模型，不依赖任何现成的机器学习库，仅使用NumPy完成矩阵运算，最终在经典的波士顿房价数据集上进行预测和评估。

1. 环境准备与数据理解

在开始编码之前，我们需要确保开发环境配置正确，并充分理解将要使用的数据集特性。

首先创建一个新的Python环境（推荐使用conda或venv），安装以下必要依赖：

pip install numpy matplotlib scikit-learn

波士顿房价数据集包含506个样本，每个样本有13个特征和1个目标值（房价中位数）。让我们先加载并探索数据：

from sklearn.datasets import load_boston
import numpy as np

boston = load_boston()
X = boston.data
y = boston.target

print(f"特征矩阵形状: {X.shape}")
print(f"目标值形状: {y.shape}")
print(f"特征名称: {boston.feature_names}")

注意：较新版本的scikit-learn已移除了波士顿房价数据集，可考虑使用fetch_california_housing()替代，或手动下载原始数据。

数据探索的几个关键点：

• 特征缩放 ：不同特征的量纲差异很大（如CRIM范围0-89，NOX范围0.3-0.9），应考虑标准化
• 异常值处理 ：部分特征如CRIM（犯罪率）存在极端值
• 特征相关性 ：某些特征如RAD和TAX高度相关，可能影响模型稳定性

2. 线性回归数学原理与实现

线性回归的核心是通过最小化预测值与真实值之间的差距，找到最优的权重参数。我们主要实现两种解法：正规方程和梯度下降。

2.1 正规方程解法

正规方程提供了一种解析解，直接通过矩阵运算得到最优参数：

θ = (XᵀX)⁻¹Xᵀy

对应的Python实现：

def normal_equation(X, y):
    # 添加偏置项（全1列）
    X_b = np.c_[np.ones((X.shape[0], 1)), X]
    
    # 计算正规方程
    theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
    
    return theta

2.2 梯度下降实现

对于大规模数据，梯度下降通常更高效。我们实现批量梯度下降：

def gradient_descent(X, y, learning_rate=0.01, n_iters=1000):
    m = X.shape[0]
    X_b = np.c_[np.ones((m, 1)), X]
    theta = np.random.randn(X_b.shape[1], 1)
    
    for _ in range(n_iters):
        gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y.reshape(-1, 1))
        theta -= learning_rate * gradients
    
    return theta

提示：实际应用中应添加早停机制（当损失变化小于阈值时停止迭代）

3. 模型训练与评估

有了算法实现后，我们需要将数据划分为训练集和测试集，训练模型并进行评估。

3.1 数据预处理

良好的数据预处理能显著提升模型性能：

from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

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

# 特征标准化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

3.2 训练模型

使用我们实现的正规方程解法训练模型：

theta = normal_equation(X_train_scaled, y_train)
print("模型参数:", theta)

3.3 评估指标实现

实现两个核心评估指标：均方误差(MSE)和决定系数(R²)

def mse(y_true, y_pred):
    return np.mean((y_true - y_pred)**2)

def r2_score(y_true, y_pred):
    ss_res = np.sum((y_true - y_pred)**2)
    ss_tot = np.sum((y_true - np.mean(y_true))**2)
    return 1 - (ss_res / ss_tot)

评估模型性能：

# 预测测试集
X_test_b = np.c_[np.ones((X_test_scaled.shape[0], 1)), X_test_scaled]
y_pred = X_test_b.dot(theta)

# 计算指标
print(f"MSE: {mse(y_test, y_pred):.2f}")
print(f"R²: {r2_score(y_test, y_pred):.2f}")

4. 模型优化与调试

基础模型建立后，我们需要考虑如何提升性能并解决常见问题。

4.1 特征工程

尝试不同的特征组合和变换：

• 多项式特征：引入特征的平方项或交叉项
• 特征选择：使用统计方法选择重要特征
• 异常值处理：识别并处理极端样本

# 示例：添加多项式特征
from sklearn.preprocessing import PolynomialFeatures

poly = PolynomialFeatures(degree=2, include_bias=False)
X_train_poly = poly.fit_transform(X_train_scaled)
X_test_poly = poly.transform(X_test_scaled)

4.2 正则化处理

为防止过拟合，可以引入L2正则化（岭回归）：

def ridge_regression(X, y, alpha=1.0):
    X_b = np.c_[np.ones((X.shape[0], 1)), X]
    m, n = X_b.shape
    I = np.eye(n)
    I[0, 0] = 0  # 不惩罚偏置项
    
    theta = np.linalg.inv(X_b.T.dot(X_b) + alpha * I).dot(X_b.T).dot(y)
    return theta

4.3 学习曲线分析

通过绘制学习曲线诊断模型问题：

import matplotlib.pyplot as plt

def plot_learning_curve(X_train, y_train, X_val, y_val):
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):
        theta = normal_equation(X_train[:m], y_train[:m])
        y_train_pred = np.c_[np.ones((m, 1)), X_train[:m]].dot(theta)
        y_val_pred = np.c_[np.ones((len(X_val), 1)), X_val].dot(theta)
        
        train_errors.append(mse(y_train[:m], y_train_pred))
        val_errors.append(mse(y_val, y_val_pred))
    
    plt.plot(train_errors, "r-+", linewidth=2, label="训练集")
    plt.plot(val_errors, "b-", linewidth=3, label="验证集")
    plt.legend()
    plt.xlabel("训练集大小")
    plt.ylabel("MSE")

5. 完整项目结构与扩展思考

一个完整的机器学习项目应该包含以下模块：

boston_housing/
├── data/                # 原始数据
├── notebooks/           # Jupyter笔记本
├── src/
│   ├── __init__.py
│   ├── data.py          # 数据加载与预处理
│   ├── models.py        # 模型实现
│   ├── evaluate.py      # 评估指标
│   └── visualize.py     # 可视化
└── README.md

进一步扩展方向：

• 模型部署 ：将训练好的模型打包为API服务
• 自动化流水线 ：使用MLflow或Airflow构建自动化训练流程
• 模型解释 ：分析各特征对预测结果的影响程度

在实际项目中，我发现特征工程往往比模型选择更能影响最终效果。比如在房价预测中，创造性地组合地理位置特征（如到市中心的距离）有时能带来意外惊喜。