最小二乘法 Python 3.7 实战：从公式推导到波士顿房价预测

1. 理解最小二乘法的数学本质

当我们面对一组散点数据时，如何找到一条最合适的直线来描述这些数据点的趋势？这正是最小二乘法要解决的核心问题。想象你是一名房地产分析师，手上有100套房子的面积和价格数据，如何用数学方法揭示这两者之间的关系？

最小二乘法的核心思想是： 最小化预测值与真实值之间的平方误差和 。用数学公式表示就是：

min Σ(y_i - ŷ_i)^2

其中y_i是真实值，ŷ_i是我们的预测值。对于一元线性回归，预测值可以表示为：

ŷ = ax + b

这里a是斜率，b是截距。我们的目标就是找到最优的a和b，使得所有数据点的预测误差平方和最小。

2. 公式推导与Python实现

2.1 数学推导过程

要最小化误差函数E(a,b) = Σ(y_i - (ax_i + b))^2，我们需要对a和b分别求偏导并令其等于0：

对b求偏导：

∂E/∂b = -2Σ(y_i - ax_i - b) = 0
=> Σy_i = aΣx_i + nb

对a求偏导：

∂E/∂a = -2Σx_i(y_i - ax_i - b) = 0
=> Σx_i y_i = aΣx_i^2 + bΣx_i

解这组方程可以得到：

a = (nΣx_i y_i - Σx_i Σy_i) / (nΣx_i^2 - (Σx_i)^2)
b = (Σy_i - aΣx_i) / n

2.2 Python代码实现

import numpy as np

class LinearRegression:
    def __init__(self):
        self.a = 0  # 斜率
        self.b = 0  # 截距
    
    def fit(self, x, y):
        n = len(x)
        sum_x = np.sum(x)
        sum_y = np.sum(y)
        sum_xy = np.sum(x * y)
        sum_x2 = np.sum(x ** 2)
        
        # 计算斜率a和截距b
        self.a = (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - sum_x ** 2)
        self.b = (sum_y - self.a * sum_x) / n
        
    def predict(self, x):
        return self.a * x + self.b

这个实现清晰地展示了最小二乘法的核心计算过程。我们通过numpy的向量化操作高效地完成了各种求和运算。

3. 波士顿房价数据集实战

3.1 数据准备与探索

波士顿房价数据集是经典的回归分析数据集，包含506个样本和13个特征。我们将使用房间数(RM)作为自变量，房价(MEDV)作为因变量。

from sklearn.datasets import load_boston
import pandas as pd
import matplotlib.pyplot as plt

# 加载数据
boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df['MEDV'] = boston.target

# 可视化RM与MEDV的关系
plt.scatter(df['RM'], df['MEDV'])
plt.xlabel('Average number of rooms per dwelling')
plt.ylabel('Median value of owner-occupied homes ($1000s)')
plt.title('Boston Housing Data')
plt.show()

从散点图可以明显看出房间数与房价之间存在正相关关系，适合用线性回归建模。

3.2 模型训练与评估

使用我们实现的最小二乘法模型：

# 准备数据
X = df['RM'].values
y = df['MEDV'].values

# 训练模型
model = LinearRegression()
model.fit(X, y)

# 预测并可视化
y_pred = model.predict(X)

plt.scatter(X, y, label='Actual')
plt.plot(X, y_pred, color='red', label='Predicted')
plt.xlabel('Number of rooms')
plt.ylabel('House price ($1000s)')
plt.legend()
plt.show()

为了评估模型性能，我们计算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)

print(f'R² score: {r2_score(y, y_pred):.3f}')

4. 与Scikit-learn的对比验证

为了验证我们实现的正确性，我们使用Scikit-learn的LinearRegression进行对比：

from sklearn.linear_model import LinearRegression as SkLinearRegression
from sklearn.model_selection import train_test_split

# 使用相同的训练测试集划分
X_train, X_test, y_train, y_test = train_test_split(
    X.reshape(-1, 1), y, test_size=0.2, random_state=42)

# 我们的实现
our_model = LinearRegression()
our_model.fit(X_train.ravel(), y_train)
our_pred = our_model.predict(X_test.ravel())
our_r2 = r2_score(y_test, our_pred)

# Scikit-learn实现
sk_model = SkLinearRegression()
sk_model.fit(X_train, y_train)
sk_pred = sk_model.predict(X_test)
sk_r2 = sk_model.score(X_test, y_test)

print(f'Our R²: {our_r2:.3f}, Sklearn R²: {sk_r2:.3f}')
print(f'Our coefficients: a={our_model.a:.3f}, b={our_model.b:.3f}')
print(f'Sklearn coefficients: a={sk_model.coef_[0]:.3f}, b={sk_model.intercept_:.3f}')

典型输出结果可能如下：

Our R²: 0.485, Sklearn R²: 0.485
Our coefficients: a=9.102, b=-34.671
Sklearn coefficients: a=9.102, b=-34.671

可以看到两者的R²分数和系数几乎完全一致，验证了我们实现的正确性。

5. 模型优化与扩展思考

虽然我们实现了一个基本的最小二乘法回归，但在实际应用中还可以考虑以下优化方向：

1. 正则化 ：加入L1/L2正则项防止过拟合
2. 多元线性回归 ：扩展到多个特征的情况
3. 梯度下降 ：对于大数据集，使用迭代优化方法
4. 特征工程 ：对特征进行标准化、多项式扩展等处理

# 多元线性回归的矩阵形式实现
class MultipleLinearRegression:
    def __init__(self):
        self.coef_ = None
        
    def fit(self, X, y):
        # 添加截距项
        X = np.column_stack([np.ones(X.shape[0]), X])
        # 最小二乘解：(X^T X)^-1 X^T y
        self.coef_ = np.linalg.inv(X.T @ X) @ X.T @ y
        
    def predict(self, X):
        X = np.column_stack([np.ones(X.shape[0]), X])
        return X @ self.coef_

这个矩阵形式的实现可以处理任意数量的特征，展示了最小二乘法更通用的表达形式。