线性回归实战-Boston 数据集回归分析

要求： 使用 scikit-learn 的 load_boston 加载波士顿房价数据集，并完成以下任务：

注：load_boston从sklearn 1.2版本之后已经被移除，所以可以采用下面的方式加载数据：

    import pandas as pd
    import numpy as np

    data_url = "http://lib.stat.cmu.edu/datasets/boston"
    raw_df = pd.read_csv(data_url, sep="\s+", skiprows=22, header=None)
    data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])
    target = raw_df.values[1::2, 2]

将数据集划分为训练集和测试集（使用 train_test_split）。
分别使用 LinearRegression、Ridge 和 Lasso 模型进行训练。
在测试集上评估每个模型的性能，使用均方误差（MSE）和 R2 分数作为指标。
比较三个模型在系数上的差异，并简要分析原因。

【1】加载数据：（在了解原数据集的基础上，我将feature_names参数中的列名，换成了中文）

# 从原始来源加载波士顿房价数据集
data_url = "http://lib.stat.cmu.edu/datasets/boston"
raw_df = pd.read_csv(data_url, sep="\s+", skiprows=22, header=None)
data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])
target = raw_df.values[1::2, 2]

# 定义中文特征名称
chinese_feature_names = [
    "犯罪率",          # CRIM: 人均犯罪率
    "住宅用地比例",    # ZN: 住宅用地比例
    "非零售商业用地比例",  # INDUS: 非零售商业用地比例
    "是否临河",        # CHAS: 是否临查尔斯河（1=是，0=否）
    "一氧化氮浓度",    # NOX: 一氧化氮浓度
    "平均房间数",      # RM: 每栋住宅的平均房间数
    "老房比例",        # AGE: 1940年前建成的房屋比例
    "距离就业中心距离",  # DIS: 到五个波士顿就业中心的加权距离
    "高速公路可达性",  # RAD: 高速公路可达性指数
    "房产税率",        # TAX: 每万美元房产的税率
    "师生比例",        # PTRATIO: 学生与教师比例
    "黑人比例",        # B: 与黑人人口比例相关的指标（有伦理争议）
    "低收入人群比例"   # LSTAT: 低收入人群比例
]

# 转换为 DataFrame，使用中文特征名称
df = pd.DataFrame(data, columns=chinese_feature_names)
df['房价'] = target  # 目标值也用中文命名

# 打印数据信息
print("\n中文特征名称:", chinese_feature_names)
print(f"数据规模: {data.shape}")
print(f"目标值规模: {target.shape}")
print("\n数据集前几行:\n")
display(df.head())

data的规模是(506, 13)，即506个样本，13个特征。

load_boston从sklearn 1.2版本之后已经被移除，官方对于波士顿数据集的加载，除了提供了上面的方式，还说明了可以采用下面2个数据集。对比整理如下：

特性	Boston Housing	California Housing	Ames Housing
样本数	506	20,640	2,930
特征数	13	8	~80（含类别变量）
目标值	房价（千美元）	房价（十万美元）	房价（美元）
伦理问题	有（特征 “B” 和研究假设争议）	无	无
数据复杂度	简单，适合初学者	中等，适合一般回归任务	复杂，需特征工程
加载方式	已移除，需从原始来源获取	fetch_california_housing()	fetch_openml(name=”house_prices”)

【2】数据集划分为训练集和测试集

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, t_test = train_test_split(data, target, test_size=0.2)

其中X_train的规模是(404, 13)， X_test的规模是(102, 13)

【3】使用 LinearRegression

# 创建线性回归模型实例
model = LinearRegression()

# 使用训练数据拟合模型
model.fit(X_train, y_train)

# 打印系数和截距
print("模型系数 (coef_):")
# 将特征名与系数对应起来
for i, name in enumerate(chinese_feature_names):
    print(f"  {name}: {model.coef_[i]:.4f}")

print(f"\n模型截距 (intercept_): {model.intercept_:.4f}")

插入一道选择题：

如果一个线性回归模型中某个特征的系数是负值，这意味着什么？

A. 该特征不重要。

B. 该特征对目标值有负面影响。

C. 该特征的单位对模型训练很重要。

D. 模型没有正确训练。

【答案】：负系数表示当该特征值增加时，目标值会减少，即两者呈负相关。

找出对房价影响最大的三个特征：

# 获取系数的绝对值，并进行排序
abs_coefs = np.abs(model.coef_)
sorted_indices = np.argsort(abs_coefs)[::-1] # 降序排列

print("对房价影响最大的三个特征:")
for i in sorted_indices[:3]:
    print(f"  {chinese_feature_names[i]}: 系数绝对值为 {abs_coefs[i]:.4f}")

【4】预测：

# 使用模型对测试集进行预测
y_pred = model.predict(X_test)

print("前10个样本的真实值 vs 预测值:")
for i in range(10):
    print(f"  真实值: {y_test[i]:.2f}, 预测值: {y_pred[i]:.2f}")

绘制预测值 y_pred 与真实值 y_test 的散点图，观察它们的分布情况:

plt.figure(figsize=(8, 6))
plt.scatter(y_test, y_pred, alpha=0.6)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', lw=2, label='完美预测线')
plt.title('预测值 vs 真实值')
plt.xlabel('真实房价')
plt.ylabel('预测房价')
plt.legend()
plt.grid(True)
plt.show()

【5】模型评估

一般评估指标就是R方和MSE

from sklearn.metrics import mean_squared_error, r2_score
# 使用model.score()计算R^2分数
r2_score_test = model.score(X_test, y_test)
print(f"测试集 R^2 分数: {r2_score_test:.4f}")

# 使用mean_squared_error()计算MSE
mse_test = mean_squared_error(y_test, y_pred)
print(f"测试集 MSE: {mse_test:.4f}")

优化idea1:特征工程 (Feature Engineering)

特征工程：创建新特征或转换现有特征来为模型提供更多信息。

对于线性回归模型，可以尝试下面几种方法：

创建交互项： 某些特征的组合可能比它们各自单独存在时对房价的影响更大。例如，高犯罪率（CRIM）与低收入人群比例（LSTAT）的乘积，可能比单独的犯罪率或低收入人群比例更能体现房价的负面影响。使用 sklearn.preprocessing.PolynomialFeatures 来自动生成多项式和交互项。
非线性特征转换： 房价与某些特征的关系可能不是线性的。例如，房价可能随平均房间数（RM）的增加呈非线性上升。可以对特征进行对数转换、平方根转换或取平方等操作，以更好地捕捉这种非线性关系。

创建多项式和交互项

78怎么来的？

from sklearn.preprocessing import PolynomialFeatures

# 创建degree=2的多项式特征
poly = PolynomialFeatures(degree=2, include_bias=False) 
X_train_poly = poly.fit_transform(X_train)
X_test_poly = poly.transform(X_test)

# 重新训练模型
model_poly = LinearRegression()
model_poly.fit(X_train_poly, y_train)
r2_score_poly = model_poly.score(X_test_poly, y_test)
print(f"使用多项式特征后的R^2分数: {r2_score_poly:.4f}")

改变参数degree会不会有更好的效果呢？下面对于参数进行了一些说明：

degree=1：这相当于一个标准的线性回归模型，没有任何多项式或交互项。它只会使用原始特征，是模型的基线。
degree=2：这是最常用的值。它会为每个特征添加平方项和所有特征的两两交互项。这通常是一个很好的起点，因为它能捕捉到大多数重要的非线性关系和特征之间的联动，同时又不会导致特征数量爆炸得太过严重。
degree=3：这个值偶尔也会使用，特别是在数据中存在更复杂的曲线关系时。但正如我们之前讨论的，它会急剧增加特征数量，因此需要更谨慎地使用，并且通常要结合正则化（如 Ridge 或 Lasso）来防止过拟合。
degree >= 4： 极少使用。当 degree 超过3时，特征数量会以指数级增长，导致模型非常容易过拟合。这会使模型在训练集上表现完美，但在新的、未见过的数据上表现极差。同时，计算成本也会变得非常高，并且模型的可解释性几乎完全丧失。

可以得出结论，degree=2就是目前最佳的效果。R2达到了0.8666

这个实验结果完美地说明了在机器学习中一个非常重要的概念：偏差-方差权衡（Bias-Variance Trade-off）。

degree=1 的模型偏差高（模型太简单，无法捕捉到数据的非线性关系），方差低（对训练数据的微小变化不敏感）。
degree=3 和 degree=4 的模型偏差低（模型足够复杂，可以完美拟合训练数据），但方差极高（对训练数据的微小变化极其敏感，导致在测试集上性能极差）。
degree=2 的模型找到了一个很好的平衡点，它的偏差和方差都处于一个相对较低的水平，使得模型既能捕捉到数据的关键模式，又具备良好的泛化能力。

优化idea2:数据预处理和标准化

虽然线性回归模型对特征缩放不敏感，但正则化模型（Ridge 和 Lasso）对特征缩放非常敏感。特征缩放可以确保所有特征在训练过程中对模型的影响力是均衡的。

标准化 (Standardization)： 使用 StandardScaler 将特征缩放到均值为0，方差为1。这是最常用的方法。

分别使用 LinearRegression、Ridge 和 Lasso 模型，但这次我们会先对数据进行标准化处理，然后对比标准化前后模型性能的变化，并深入分析。

步骤1:标准化数据

from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.metrics import mean_squared_error, r2_score
import numpy as np

# 假设 X_train, X_test, y_train, y_test 已经通过 train_test_split 划分好

# 实例化 StandardScaler
scaler = StandardScaler()

# 在训练集上拟合 (fit) 并转换 (transform)
X_train_scaled = scaler.fit_transform(X_train)

# 使用同一个 scaler 转换测试集
X_test_scaled = scaler.transform(X_test)

步骤2:使用标准化后的数据 X_train_scaled 和 X_test_scaled 分别训练 LinearRegression、Ridge 和 Lasso 模型，并评估它们的性能。

# 初始化模型
models = {
    "LinearRegression": LinearRegression(),
    "Ridge": Ridge(alpha=1.0),
    "Lasso": Lasso(alpha=0.1)
}

# 存储结果
results = {}

print("--- 标准化后的模型性能评估 ---")
for name, model in models.items():
    # 训练模型
    model.fit(X_train_scaled, y_train)

    # 预测
    y_pred = model.predict(X_test_scaled)

    # 评估
    mse = mean_squared_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    
    results[name] = {"MSE": mse, "R2 Score": r2}
    
    print(f"\n模型: {name}")
    print(f"  测试集 MSE: {mse:.4f}")
    print(f"  测试集 R^2 分数: {r2:.4f}")

print("\n--- 性能对比（标准化前 vs 标准化后）---")
# 你可以把之前未标准化的 LinearRegression 结果也放进来对比
# 假设之前的 LinearRegression 结果是 r2_score_unscaled 和 mse_unscaled
# print(f"LinearRegression (未标准化): R^2 = {r2_score_unscaled:.4f}, MSE = {mse_unscaled:.4f}")
for name, metrics in results.items():
    print(f"{name} (标准化后): R^2 = {metrics['R2 Score']:.4f}, MSE = {metrics['MSE']:.4f}")

可以看到效果一般。

优化方案：多项式特征 + 标准化 + 正则化

# 1. 创建多项式特征 (degree=2)
from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures(degree=2, include_bias=False)
X_train_poly = poly.fit_transform(X_train)
X_test_poly = poly.transform(X_test)

# 2. 标准化新特征
from sklearn.preprocessing import StandardScaler
scaler_poly = StandardScaler()
X_train_poly_scaled = scaler_poly.fit_transform(X_train_poly)
X_test_poly_scaled = scaler_poly.transform(X_test_poly)

# 3. 使用正则化模型进行训练和评估
# 我们可以尝试不同的 alpha 值，这里先用一个默认值
ridge_poly_model = Ridge(alpha=1.0)
ridge_poly_model.fit(X_train_poly_scaled, y_train)

y_pred_ridge_poly = ridge_poly_model.predict(X_test_poly_scaled)
r2_ridge_poly = r2_score(y_test, y_pred_ridge_poly)
mse_ridge_poly = mean_squared_error(y_test, y_pred_ridge_poly)

print(f"多项式特征 + 标准化 + Ridge 模型性能:")
print(f"  测试集 R^2 分数: {r2_ridge_poly:.4f}")
print(f"  测试集 MSE: {mse_ridge_poly:.4f}")

这非常棒！R2 值从 0.8666 提升到了 0.8740，这说明在多项式特征的基础上引入标准化和正则化是有效的。尽管提升看起来不大。

多项式特征提升了模型的表达能力，而正则化则防止了模型在复杂化后过拟合。

优化 Idea 3: 超参数调优 (Hyperparameter Tuning)

Ridge 和 Lasso 模型中使用的 alpha 值（例如 Ridge(alpha=1.0)）是手动设定的，这可能不是最优解。alpha 控制着正则化的强度，不同的 alpha 值会显著影响模型的表现。

Scikit-learn 提供了 RidgeCV 和 LassoCV，它们是内置了交叉验证的 Ridge 和 Lasso 模型。你只需要提供一个 alpha 值的列表，模型就会自动为你找到最佳的 alpha，并用它来训练模型。

from sklearn.linear_model import RidgeCV, LassoCV
from sklearn.model_selection import KFold
import numpy as np

# 假设 X_train_poly_scaled 和 y_train 已经准备就绪
# 定义一组可能的 alpha 值
alphas = np.logspace(-3, 2, 50) 

# 定义交叉验证策略
cv = KFold(n_splits=5, shuffle=True, random_state=42)

print("--- 增加迭代次数解决 ConvergenceWarning ---")

# 使用 RidgeCV，并增加 max_iter
# 注意：在 RidgeCV 中，max_iter 参数是直接传递给 fit 方法的
# 或者在实例化时通过其他方式，但最简单的是使用一个常规的 Ridge 模型与 GridSearchCV
# 让我们换成一个更通用的方法，使用 GridSearchCV 来处理超参数调优
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import Ridge, Lasso

# 为 Ridge 模型创建一个参数网格
param_grid_ridge = {'alpha': alphas}

# 使用 GridSearchCV 进行交叉验证
grid_search_ridge = GridSearchCV(Ridge(max_iter=100000), param_grid_ridge, cv=cv)
grid_search_ridge.fit(X_train_poly_scaled, y_train)

# 打印最佳 alpha 值和性能
print(f"GridSearchCV for Ridge 找到的最佳 alpha 值: {grid_search_ridge.best_params_['alpha']:.4f}")
y_pred_ridge_grid = grid_search_ridge.best_estimator_.predict(X_test_poly_scaled)
r2_ridge_grid = r2_score(y_test, y_pred_ridge_grid)
print(f"GridSearchCV for Ridge 最佳 alpha 下的 R^2 分数: {r2_ridge_grid:.4f}")

# 使用 LassoCV，max_iter 是一个直接支持的参数
# 让我们保持 LassoCV，因为它在这里工作得很好
lasso_cv_fixed = LassoCV(alphas=alphas, cv=cv, max_iter=100000)
lasso_cv_fixed.fit(X_train_poly_scaled, y_train)

# 打印最佳 alpha 值和性能
print(f"\nLassoCV 找到的最佳 alpha 值: {lasso_cv_fixed.alpha_:.4f}")
y_pred_lasso_cv_fixed = lasso_cv_fixed.predict(X_test_poly_scaled)
r2_lasso_cv_fixed = r2_score(y_test, y_pred_lasso_cv_fixed)
print(f"LassoCV 最佳 alpha 下的 R^2 分数: {r2_lasso_cv_fixed:.4f}")

目前最佳是0.8770

可视化特征

线性回归及其变种（Ridge、Lasso）对异常值非常敏感。一个或几个离群点可能就会显著影响模型的系数，从而导致预测偏差。

绘制每个特征与房价之间的散点图。特别是，关注那些看起来与主要数据分布不一致的点。

import seaborn as sns
# 获取特征列名（排除目标列 '房价'）
feature_columns = [col for col in df.columns if col != '房价']

# 设置图表样式
sns.set_style("whitegrid")
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False     # 用来正常显示负号

# 遍历所有特征并绘制散点图
fig, axes = plt.subplots(nrows=7, ncols=2, figsize=(15, 30))
fig.suptitle('各特征与房价的散点图', fontsize=16, y=0.92)

# 将 axes 展平，方便迭代
axes = axes.flatten()

for i, col in enumerate(feature_columns):
    sns.scatterplot(x=df[col], y=df['房价'], ax=axes[i], alpha=0.6)
    axes[i].set_title(f'第{i}个特征:{col} vs 房价')
    axes[i].set_xlabel(col)
    axes[i].set_ylabel('房价')

# 隐藏最后一个空白的子图
axes[13].axis('off')

plt.tight_layout(rect=[0, 0, 1, 0.9])
plt.show()

犯罪率和房价存在明显的负相关。随着犯罪率的增加，房价总体呈下降趋势。图中的一些异常点，如高犯罪率但房价不低的样本，可能值得关注。
住宅用地比例与房价关系不太明确，但略呈正相关。当住宅用地比例较高时，房价似乎有更高的分布，尤其是在 ZN 接近 100% 的地方，房价普遍较高。
非零售商业用地比例与房价存在明显的负相关。商业用地比例越高，房价越低。
是否临河是一个二元变量（0或1）。临河（CHAS=1）的房屋，房价的分布总体上高于不临河（CHAS=0）的房屋。
一氧化氮浓度与房价存在明显的负相关。空气污染（一氧化氮浓度）越高，房价越低。
平均房间数与房价存在非常强的正相关。房间数越多，房价越高。这几乎是数据集中最强的正相关关系之一。
老房比例与房价存在一定的负相关。房屋越旧，房价越可能偏低。但也存在一些高房价的老房子，可能说明其位于好地段或经过翻修。
距离就业中心距离与房价存在正相关。与就业中心的距离越大，房价反而越高。这可能是一个反直觉的结论，但它可能表明，波士顿的房价在郊区（距离就业中心远）反而更高，因为这些地方可能拥有更好的环境和学区，吸引高收入人群。
高速公路可达性与房价关系不明显，但RAD值较小的区域，房价分布更广，既有高价房也有低价房。而RAD值较高的区域（如24），房价普遍较低。高速公路可达性指数较高的区域可能意味着更嘈杂或工业化的环境，从而对房价有负面影响。
房产税率与房价存在明显的负相关。税率越高，房价越低。图中的垂直散点带表明，在某个特定的税率下，房价的分布非常广，这可能是数据本身的特性。
师生比例与房价存在负相关。师生比例越高，房价越低。这通常与学区质量相关，较低的师生比例（即更好的教育资源）是房价的正面因素。图中有一些离群点，例如师生比例很高但房价也高的点，值得进一步研究。
黑人比例与房价存在正相关。黑人比例相关的指标越高，房价越高。这个特征是数据集中最具争议性的特征，原始数据集的背景是种族隔离和区域贫富差距。图中的正相关关系可能反映了当时社会经济结构中某些特定群体的集中和相应的房价水平。在现代数据分析中，这个特征通常会被谨慎处理或移除。
低收入人群比例与房价存在非常强的负相关。LSTAT 越高，房价越低。低收入人群比例是影响房价的最重要负面因素之一。图中的非线性关系非常明显，它呈现出一条下凹的曲线，说明LSTAT在较低水平时，房价对它的变化不那么敏感，但在LSTAT值较高时，房价会急剧下降。

机器学习算法2-线性回归-波士顿房价问题