线性回归实战：从MSE推导到Python代码实现（附完整数据集）

在机器学习领域，线性回归就像"Hello World"一样经典而重要。但很多初学者都会遇到这样的困境：看懂了数学公式，却不知道如何用代码实现；理解了理论概念，却无法将其应用到真实数据集上。本文将彻底解决这个痛点，带你从MSE的数学推导开始，一步步实现完整的线性回归模型，最后用Python代码在真实房价数据集上验证理论。

1. 线性回归与MSE的数学本质

线性回归的核心思想非常简单：找到一条直线（或超平面），使得所有数据点到这条直线的垂直距离之和最小。这个"距离"在数学上通过**均方误差（MSE）**来量化：

$$ MSE = \frac{1}{n}\sum_{i=1}^n(y_i - \hat{y_i})^2 $$

为什么选择平方而不是绝对值？这背后有深刻的统计学原理：

1. 最大似然估计视角：假设误差服从正态分布，最大化似然函数等价于最小化MSE
2. 数学便利性：平方函数处处可导，便于优化求解
3. 对大误差的惩罚：平方放大了大误差的影响，使模型更关注严重错误

在实际房价预测中，MSE的单位是"万元平方"，虽然解释性不如绝对误差，但作为优化目标非常有效

2. 从零推导MSE损失函数

让我们从概率角度完整推导MSE的由来：

2.1 基本假设

• 线性关系：$y = X\theta + \epsilon$
• 误差项$\epsilon$服从$N(0,\sigma^2)$，即均值为0的正态分布

2.2 似然函数构建

单个样本的概率密度：

$$ p(\epsilon_i) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\epsilon_i^2}{2\sigma^2}\right) $$

由于误差独立同分布，整体似然函数为：

$$ L(\theta) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y_i - x_i\theta)^2}{2\sigma^2}\right) $$

2.3 对数似然与简化

取对数后得到：

$$ \ln L(\theta) = -\frac{n}{2}\ln(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n(y_i - x_i\theta)^2 $$

最大化似然等价于最小化：

$$ J(\theta) = \sum_{i=1}^n(y_i - x_i\theta)^2 $$

这就是MSE的原始形式！

3. Python实现：纯NumPy版本

理解了数学原理后，我们先用NumPy实现最基础的版本：

import numpy as np

class LinearRegression:
    def __init__(self):
        self.weights = None
    
    def fit(self, X, y):
        # 添加偏置项
        X = np.c_[np.ones(X.shape[0]), X]
        
        # 解析解：(X^T X)^-1 X^T y
        self.weights = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
    
    def predict(self, X):
        X = np.c_[np.ones(X.shape[0]), X]
        return X.dot(self.weights)

关键点说明：

• np.c_用于添加全1列作为偏置项
• 直接使用解析解公式，避免迭代计算
• 矩阵运算比循环效率高几个数量级

4. 在房价数据集上的实战

我们使用波士顿房价数据集进行验证：

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split

# 加载数据
boston = load_boston()
X, y = boston.data, boston.target

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

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

# 评估
train_pred = model.predict(X_train)
test_pred = model.predict(X_test)

print(f'Train MSE: {np.mean((y_train - train_pred)**2):.2f}')
print(f'Test MSE: {np.mean((y_test - test_pred)**2):.2f}')

典型输出结果：

Train MSE: 21.64
Test MSE: 24.29

5. 与Scikit-learn的实现对比

为了验证我们的实现，与行业标准库进行对比：

指标	我们的实现	Scikit-learn
训练MSE	21.64	21.64
测试MSE	24.29	24.29
运行时间(ms)	2.1	1.8

实现对比代码：

from sklearn.linear_model import LinearRegression as SkLinearRegression

sk_model = SkLinearRegression()
sk_model.fit(X_train, y_train)

sk_train_pred = sk_model.predict(X_train)
sk_test_pred = sk_model.predict(X_test)

print('\nScikit-learn结果:')
print(f'Train MSE: {np.mean((y_train - sk_train_pred)**2):.2f}')
print(f'Test MSE: {np.mean((y_test - sk_test_pred)**2):.2f}')

6. 进阶技巧与优化

6.1 数值稳定性优化

解析解中的矩阵求逆可能不稳定，可以改用SVD分解：

def fit_svd(self, X, y):
    X = np.c_[np.ones(X.shape[0]), X]
    U, s, Vt = np.linalg.svd(X, full_matrices=False)
    self.weights = Vt.T @ np.diag(1/s) @ U.T @ y

6.2 特征缩放的重要性

对于不同量纲的特征，应先进行标准化：

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

model.fit(X_train_scaled, y_train)

6.3 梯度下降实现

当特征维度很高时，解析解计算成本高，可用梯度下降：

def fit_gd(self, X, y, lr=0.01, epochs=1000):
    X = np.c_[np.ones(X.shape[0]), X]
    n_samples, n_features = X.shape
    self.weights = np.zeros(n_features)
    
    for _ in range(epochs):
        grad = (2/n_samples) * X.T @ (X @ self.weights - y)
        self.weights -= lr * grad

7. 项目实践建议

在实际房价预测项目中，还需要考虑：

1. 特征工程：

   • 处理缺失值
   • 异常值检测
   • 特征交叉

2. 模型诊断：

   # 残差分析
   residuals = y_test - test_pred
   plt.scatter(test_pred, residuals)
   plt.axhline(y=0, color='r', linestyle='-')
   

3. 正则化：

   • 岭回归(L2正则)
   • Lasso回归(L1正则)

完整项目代码和数据集已上传至GitHub仓库，包含：

• 数据清洗脚本
• 特征工程示例
• 多种实现方式对比
• Jupyter Notebook教程