从零实现线性回归：用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 模型假设验证

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

1. 线性关系：特征和目标之间确实存在线性关系
2. 同方差性：残差的方差应该是常数
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}")