保姆级教程：用Python手撸线性回归，从波士顿房价预测到模型评估（附完整代码）-编程实验室

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

在机器学习领域，线性回归堪称算法界的"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

进一步扩展方向：