当前位置：首页 > news >正文

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

news 2026/5/31 11:35:15

从零实现线性回归：用Python揭秘房价预测的数学之美

在机器学习的世界里，线性回归就像是一把打开预测大门的金钥匙。它不仅是最基础的算法之一，更是理解更复杂模型的必经之路。想象一下，你能够通过几个简单的数学公式，就能预测一座房子的价格——这正是线性回归的魅力所在。本文将带你从数学原理出发，手把手实现一个完整的线性回归模型，并用它来预测波士顿房价。不同于简单的调用sklearn，我们会深入算法内部，用NumPy一步步构建模型，让你真正理解机器学习的"黑箱"里发生了什么。

1. 线性回归的数学基础

线性回归的核心思想非常简单：找到一条直线（在更高维度是超平面），使得这条直线能够最好地"拟合"数据点。用数学语言来说，就是找到一组参数θ，使得预测值ŷ与真实值y之间的差距最小。

1.1 模型公式解析

多元线性回归的预测公式可以表示为：

ŷ = θ₀ + θ₁x₁ + θ₂x₂ + ... + θₙxₙ

其中：

ŷ是预测值
θ₀是偏置项（截距）
θ₁到θₙ是各个特征的权重系数
x₁到xₙ是输入特征

为了简化计算，我们通常会添加一个全为1的特征x₀，这样公式可以写成向量形式：

ŷ = hθ(x) = θᵀx

1.2 损失函数与优化目标

为了衡量模型的好坏，我们需要定义一个损失函数。在线性回归中，最常用的是均方误差(MSE)：

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

我们的目标就是找到一组θ，使得MSE最小化。这个优化问题有解析解，称为"正规方程"：

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

其中：

X是特征矩阵（包含x₀=1的列）
y是目标值向量

2. 数据准备与预处理

2.1 加载波士顿房价数据集

波士顿房价数据集是机器学习中的经典数据集，包含506个样本，每个样本有13个特征和1个目标值（房价中位数）。

from sklearn.datasets import load_boston boston = load_boston() X = boston.data y = boston.target

2.2 数据标准化

为了确保不同尺度的特征对模型有相同的影响，我们需要对数据进行标准化处理：

from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_scaled = scaler.fit_transform(X)

2.3 添加偏置项

记得我们在数学推导中添加的x₀=1吗？现在需要在代码中实现：

# 添加全1列作为偏置项 X_b = np.c_[np.ones((len(X_scaled), 1)), X_scaled]

3. 实现线性回归

3.1 正规方程解法

根据前面的数学推导，我们可以直接实现正规方程：

def linear_regression(X, y): # 计算θ = (XᵀX)⁻¹Xᵀy theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y) return theta

3.2 处理矩阵不可逆情况

在实际应用中，XᵀX可能不是满秩矩阵，导致无法求逆。我们可以使用伪逆来避免这个问题：

theta = np.linalg.pinv(X.T.dot(X)).dot(X.T).dot(y)

3.3 完整实现

将上述步骤整合成一个完整的类：

class LinearRegression: def __init__(self): self.theta = None def fit(self, X, y): X_b = np.c_[np.ones((len(X), 1)), X] self.theta = np.linalg.pinv(X_b.T.dot(X_b)).dot(X_b.T).dot(y) def predict(self, X): X_b = np.c_[np.ones((len(X), 1)), X] return X_b.dot(self.theta)

4. 模型评估与优化

4.1 常用评估指标

除了MSE，我们还可以计算R²分数来评估模型性能：

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)

4.2 交叉验证

为了更可靠地评估模型，我们应该使用交叉验证：

from sklearn.model_selection import cross_val_score lr = LinearRegression() mse_scores = -cross_val_score(lr, X_scaled, y, scoring='neg_mean_squared_error', cv=5) r2_scores = cross_val_score(lr, X_scaled, y, scoring='r2', cv=5)

4.3 特征重要性分析

线性回归的一个优势是我们可以直接查看各个特征的权重：

feature_importance = pd.DataFrame({ 'Feature': ['Bias'] + list(boston.feature_names), 'Weight': lr.theta }).sort_values('Weight', key=abs, ascending=False)

5. 实际应用中的注意事项

5.1 特征相关性检查

在应用线性回归前，应该检查特征之间的相关性。高度相关的特征会导致矩阵接近奇异，影响模型稳定性。

corr_matrix = pd.DataFrame(X_scaled).corr()

5.2 异常值处理

线性回归对异常值敏感，可以考虑使用以下方法：

可视化检测（箱线图、散点图）
使用RobustScaler代替StandardScaler
考虑岭回归或Lasso回归

5.3 模型假设验证

线性回归有几个重要假设需要验证：

线性关系：特征和目标之间确实存在线性关系
同方差性：残差的方差应该是常数
正态性：残差应该近似正态分布

可以通过残差图来验证这些假设：

y_pred = lr.predict(X_scaled) residuals = y - y_pred plt.scatter(y_pred, residuals) plt.axhline(y=0, color='r', linestyle='-') plt.xlabel('Predicted values') plt.ylabel('Residuals')

6. 扩展与进阶

6.1 梯度下降实现

除了正规方程，我们还可以用梯度下降来求解线性回归：

def gradient_descent(X, y, learning_rate=0.01, n_iters=1000): m = len(y) theta = np.zeros(X.shape[1]) for _ in range(n_iters): gradients = 2/m * X.T.dot(X.dot(theta) - y) theta -= learning_rate * gradients return theta

6.2 正则化方法

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

# 岭回归 theta_ridge = np.linalg.inv(X.T.dot(X) + alpha*np.eye(X.shape[1])).dot(X.T).dot(y)

6.3 多项式回归

当数据不是线性关系时，可以通过添加多项式特征来增强模型表现：

from sklearn.preprocessing import PolynomialFeatures poly = PolynomialFeatures(degree=2) X_poly = poly.fit_transform(X_scaled)

7. 完整代码示例

以下是整合了所有功能的完整实现：

import numpy as np from sklearn.datasets import load_boston from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split # 加载数据 boston = load_boston() X, y = boston.data, boston.target # 数据预处理 scaler = StandardScaler() X_scaled = scaler.fit_transform(X) X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2) # 线性回归实现 class LinearRegression: def __init__(self): self.theta = None def fit(self, X, y): X_b = np.c_[np.ones((len(X), 1)), X] self.theta = np.linalg.pinv(X_b.T.dot(X_b)).dot(X_b.T).dot(y) def predict(self, X): X_b = np.c_[np.ones((len(X), 1)), X] return X_b.dot(self.theta) def mse(self, X, y): y_pred = self.predict(X) return np.mean((y - y_pred)**2) def r2(self, X, y): y_pred = self.predict(X) ss_res = np.sum((y - y_pred)**2) ss_tot = np.sum((y - np.mean(y))**2) return 1 - (ss_res / ss_tot) # 训练和评估模型 lr = LinearRegression() lr.fit(X_train, y_train) print(f"MSE on test set: {lr.mse(X_test, y_test):.2f}") print(f"R² on test set: {lr.r2(X_test, y_test):.2f}")

查看全文

http://www.zskr.cn/news/1433802.html