线性回归+梯度下降实现PM2.5预测：机器学习大作业源码解析-开发者社区

简介：一份用于机器学习课程大作业的Python项目，收集合肥地区过去一年的月平均空气质量数据，以PM2.5为预测目标，构建线性回归模型并预测后续某月的数值。项目采用矩阵形式的线性回归和梯度下降法求解参数，完整实现了从数据读取、特征拼接、模型训练到结果评估的闭环流程，适合正在完成机器学习作业或入门回归算法的高校学生参考。压缩包共21个文件，主要由12个CSV数据文件（含训练集、测试集、特征拼接表及预测结果）、3个Python源码文件、1个模型参数npy文件以及说明文档、示例图片组成，整体大小约2.58MB，目录结构清晰，下载后即可快速运行。目前已有280人学习使用，通过该资源可以掌握线性回归的矩阵运算实现思路，并可直接复用其数据预处理和梯度下降代码开展其他空气质量预测实验。

1. 机器学习大作业：基于线性回归的 PM2.5 预测源码到底能做什么

如果你是正在找机器学习课程设计参考的在校生，或者刚接触回归任务、想看看一份不调 sklearn 的线性回归怎么写，这份「基于线性回归的 PM2.5 预测」源码包值得花半小时拆一遍。它的核心不是调包预测，而是用纯 Python 实现矩阵形式的线性回归和梯度下降——从读取合肥地区历史空气质量数据开始，构造特征矩阵，训练模型，再用训练好的参数预测未来某个月的 PM2.5 值。包里既有完整源码和训练/测试 CSV，也有训练好的 model.npy 参数和可视化图片，能帮你从数据格式到梯度更新公式完整走一遍，而不是只看到一个 predict 函数。下面按我实际拆包阅读的顺序，把数据构造、模型训练、输出评估和踩坑点一一讲清楚。

2. 先看清数据链路：train.csv 与五份中间矩阵的关系

拿到一个机器学习项目源码，我习惯先不看算法，把 CSV 文件之间的关系理清楚。因为线性回归本身不难，难的是数据怎么从原始表变成模型能吃进去的 X 和 y。这份包里 CSV 文件很多，初看有点乱，但把它们的角色梳理出来后，整套流程就很清晰了。

2.1 原始数据与派生数据：谁是输入，谁是中间产物

先区分两类文件：一类是原始输入，如 train.csv、test.csv、testdata.csv；另一类是模型训练过程中生成的中间矩阵，如 arrayx.csv、concatenateX.csv、arrayy.csv、x_t.csv、listx.csv、predict.csv、ans.csv。train.csv 是合肥地区过去一段时间的空气质量记录，每一行对应一个时间点的观测值，字段一般包含 PM2.5 浓度、日期或月份编号等。test.csv 是用于预测的输入，比如今年的某几个月；testdata.csv 可能是 test.csv 的另一种格式版本，作用等价。

接下来看中间矩阵。机器学习里，线性回归要求输入是二维矩阵 X，每一行是一个样本，每一列是一个特征；输出是一维向量 y，每一行对应样本的真实值。arrayx.csv 就是构造出来的原始特征矩阵，假设我们把过去 12 个月作为特征、当前月作为标签，那 arrayx.csv 的每一行就是一组历史月份数据。concatenateX.csv 看名字就知道是在做纵向拼接——把多段窗口的特征拼成完整矩阵。这样做的原因很直接：训练样本不够时，单个月份无法构成足够多的特征，滑窗切分可以扩大样本量。

import pandas as pd import numpy as np train = pd.read_csv('train.csv') print(train.shape) # 例如 (365, 5)，365 行按天记录 print(train.head(3)) # 确认列名和时间跨度

这里 train.shape 能快速确认数据量。如果原始表是月度平均，一年只有 12 行，直接做回归样本量太少，模型泛化能力会很差。我一般会先看一眼前三行，搞清楚列名是中文还是英文、日期格式是否统一。

# 滑窗构造特征矩阵：用过去 window 个月预测下一个月 def create_dataset(data, window=12): X, y = [], [] for i in range(len(data) - window): X.append(data[i:i+window]) # 过去 window 个值作为特征 y.append(data[i+window]) # 下一个值作为标签 return np.array(X), np.array(y) arrayx, arrayy = create_dataset(train['pm25'].values, window=12) print(arrayx.shape) # (n-12, 12)，每行 12 个连续历史值 print(arrayy.shape) # (n-12,)，对应第 13 个月的观测值

这段代码的逻辑是：把按时间排列的 PM2.5 序列切成等长窗口，前 12 个值作为特征，第 13 个值作为标签，然后窗口整体后移一位。这样做的结果是样本数从 n 变成了 n-12，虽然每个样本的特征维度固定为 12，但样本量显著增加，而且保留了时间上的先后依赖。window 是核心超参数，选 12 通常对应「用过去一年预测下一个月」的课程语义，你也可以按实际数据长度调成 6 或 24。

2.2 训练集与测试集的特征对齐

predict.csv 是模型对测试数据的预测结果，ans.csv 则可能是最终答案或预测值的后处理版本。x_t.csv 看命名是测试集的特征矩阵（X_test），即把 test.csv 按同样的滑窗逻辑处理后得到的新矩阵。listx.csv 我理解是用来记录特征对应的时间标签列表，比如第几行对应哪个月份，方便把预测值映射回真实时间轴。

xt = pd.read_csv('x_t.csv') print(xt.shape) # 行数是测试样本数，列数必须和训练特征维度一致

特征对齐是最容易翻车的地方：训练特征维度是 12，测试特征也必须是 12。如果训练用过去 12 个月，测试却想预测今年第一个月，你会发现根本没有前 12 个月的数据——这就是为什么课程大作业通常预测「中间某个月」或者需要先回溯补足历史。

s_gra.csv 和 sampleSubmission.csv 属于辅助文件。sampleSubmission.csv 是课程样例提交格式，告诉你怎么组织答案；s_gra.csv 可能是用于绘制趋势对比图的数据（smooth + graph 的缩写），把真实值和预测值放在同一时间轴上做可视化验证。image.png 和 demo.jpg 就是运行后生成的曲线图，能直观看出模型拟合效果。

到这里，数据链路的全貌就清楚了：train.csv → 滑窗构造 arrayx/arrayy → 拼接得到 concatenateX → 训练完成后用同样的窗口处理 test.csv 得到 x_t.csv → 模型输出 predict.csv → 与真实值对比得到 ans.csv 和 s_gra.csv。这个结构的优点是每一步中间结果都可审计，适合写进大作业报告；缺点是手工生成的中间文件多，容易混淆。

3. 梯度下降落地：矩阵形式线性回归的实现与参数选择

数据准备好之后，核心就是训练模型。这份源码选的是「线性回归模型 + 矩阵模型 + 梯度下降公式」，也就是不依赖 sklearn 的 LinearRegression，而是自己写出参数更新过程。这一步既是大作业的加分点，也是理解机器学习本质的关键。

3.1 损失函数与梯度更新的矩阵写法

线性回归的假设函数是 h(X) = Xθ，其中 X 是 m×n 的矩阵，θ 是 n×1 的参数向量。为了让截距项也能学出来，通常会在 X 的最左边加一列全 1，这一列就是 bias 项。损失函数用均方误差 MSE，写成矩阵形式是：

def compute_loss(X, y, theta): m = len(y) pred = X.dot(theta) loss = 1 / (2 * m) * np.sum((pred - y) ** 2) return loss

这里除以 2m 是为了后续求导时消掉系数 2，属于课程里常用的约定写法。梯度下降的更新公式是 θ = θ - α * (1/m) * X^T(Xθ - y)，其中 α 是学习率。用矩阵乘法一次算出所有样本的梯度，比逐样本更新快得多。

def gradient_descent(X, y, theta, alpha, epochs): m = len(y) for i in range(epochs): gradient = X.T.dot(X.dot(theta) - y) / m theta = theta - alpha * gradient if i % 100 == 0: loss = compute_loss(X, y, theta) print(f"epoch {i}, loss = {loss:.4f}") return theta

这段代码的核心是 gradient 那一行：X.T.dot(X.dot(theta) - y) 是 X^T 乘以残差向量，除以 m 得到平均梯度。单次迭代的运算量是 O(mn)，矩阵乘法由 numpy 底层加速，对这份数据量完全够用。epochs 和 alpha 是主要调节参数。

3.2 学习率与迭代次数怎么定

alpha 和 epochs 是梯度下降最敏感的两个超参数，我拆项目时会优先看这两个值有没有在源码里写死。如果写死了，先把代码跑一遍，观察 loss 的打印输出再决定是否调整。

alpha = 0.01 # 学习率，太大容易震荡，太小收敛慢 epochs = 1000 # 迭代次数，看 loss 是否已经平稳 theta = np.zeros((concatenateX.shape[1], 1)) # 参数初始化为 0 向量 theta = gradient_descent(concatenateX, arrayy.values.reshape(-1, 1), theta, alpha, epochs)

学习率的选择直接影响训练是否收敛。alpha=0.01 在特征值范围较小（比如 PM2.5 浓度在 0~300 之间）时通常能稳定收敛；如果 loss 前几百次迭代就在震荡甚至变大，说明 alpha 过大，要降到 0.003 或 0.001。反过来，如果 loss 下降得非常慢，1000 次迭代后还在明显递减，那就调大 alpha 或增加 epochs。判断收敛的标准不是迭代次数本身，而是 loss 曲线是否进入平台期。

np.save('model.npy', theta) # 保存训练好的参数

model.npy 就是训练产物，之后做预测时直接加载它，不需要重新训练。如果换了数据集，参数大小（维度）必须和新的特征矩阵匹配，否则加载后做矩阵乘法会直接报错。

3.3 为什么用矩阵形式而不是调用 sklearn

课程大作业里用原生梯度下降而不是 sklearn 的 LinearRegression，和成绩评价标准直接相关。老师一般会看重你是否理解模型内部机制：矩阵形式能把损失函数、梯度公式、参数更新三条核心链条清晰地展示出来。用 sklearn 三行代码解决问题，报告里能写的东西就少很多，答辩时也容易在梯度推导上被追问。

另一个原因是 sklearn 的 LinearRegression 默认使用最小二乘法（解析解），不是梯度下降，而这门课的语义明确要求梯度下降公式。解析解在特征维度不高时更快更准，但写不出「参数如何一步步更新」的过程。当然，如果你只是为了快速对比结果，可以在验证阶段用 sklearn 算一遍作参考，但主模型最好保留手写梯度下降的版本。

# 参考：sklearn 的解析解，用于和梯度下降结果做对比 from sklearn.linear_model import LinearRegression lr = LinearRegression() lr.fit(concatenateX[:, 1:], arrayy) # 注意 sklearn 自动处理截距 print(lr.coef_, lr.intercept_)

不过这份资源既然以「不调包实现」为核心，我从实际使用出发，建议你以手写版本为主、sklearn 只作为验证手段——我在拆项目时通常就是这么干的。

4. 从 model.npy 到 predict.csv：预测、评估与可视化验证

训练出模型只是第一步，能对测试集做预测并评估效果，才算完整闭环。这一章把预测、评估脚本和可视化三者串起来讲，同时给出参数怎么改的建议。

4.1 加载模型参数并预测

有了 model.npy，预测阶段的核心是「复刻训练时的特征变换」。已经处理好的 x_t.csv 就是测试特征矩阵，直接用训练好的 theta 做矩阵乘法即可。

import numpy as np theta = np.load('model.npy') x_t = pd.read_csv('x_t.csv').values # 保证测试特征和训练特征列数一致 assert x_t.shape[1] == theta.shape[0], "特征维度不匹配，请检查模型参数" # 预测 = 特征矩阵 × 参数向量 predictions = x_t.dot(theta) np.savetxt('predict.csv', predictions, delimiter=',', fmt='%.2f') print(predictions[:5])

这里最关键的是 assert 那行：训练时构造的 concatenateX 有多少列（包括 bias 列），theta 就应有几行，测试特征矩阵必须完全一致。如果预测时报 shape mismatch，十有八九是预测阶段忘了加 bias 列，或者测试集用了不同的滑窗长度。predict.csv 保存时用 fmt='%.2f' 限制小数位，能让输出文件更干净，同时避免浮点噪声干扰提交效果。

4.2 评估脚本：R² 与 RMSE 怎么算

evalu.py 这个脚本负责计算预测值和真实值之间的差距。评价指标一般有两个：均方根误差 RMSE 和决定系数 R²。RMSE 反映平均误差的绝对值，R² 反映模型解释了多少方差。课程报告里通常两个都要写。

# evalu.py 的核心逻辑 y_true = pd.read_csv('ans.csv').values.flatten() y_pred = pd.read_csv('predict.csv').values.flatten() # RMSE rmse = np.sqrt(np.mean((y_true - y_pred) ** 2)) # R² ss_res = np.sum((y_true - y_pred) ** 2) ss_tot = np.sum((y_true - np.mean(y_true)) ** 2) r2 = 1 - ss_res / ss_tot print(f"RMSE = {rmse:.4f}") print(f"R² = {r2:.4f}")

RMSE 的单位和 PM2.5 浓度一致，比如 15.23 表示平均偏差约 15μg/m³；R² 越接近 1 越好，0.7 以上在月度均值预测里已经算不错。注意 ans.csv 必须是和 predict.csv 等长的真实值序列，对齐错位会导致评估结果完全失真——这个问题我在第 5 章会细说。

4.3 可视化验证：s_gra.csv 与两张图片

单独看数字不够直观，课程报告里放一张真实值和预测值的对比折线图，说服力会强很多。s_gra.csv 应该就是为画图准备的：一列是时间序号，一列是真实值，一列是预测值。image.png 和 demo.jpg 是已生成的示例图。

import matplotlib.pyplot as plt sgra = pd.read_csv('s_gra.csv') plt.figure(figsize=(12, 5)) plt.plot(sgra['time'], sgra['true'], label='True', marker='o') plt.plot(sgra['time'], sgra['pred'], label='Pred', marker='x') plt.xlabel('Month') plt.ylabel('PM2.5') plt.title('PM2.5 Prediction vs Ground Truth') plt.legend() plt.grid(True) plt.savefig('comparison.png', dpi=150)

画图的核心目的是看趋势是否一致。如果预测曲线整体滞后真实值一个相位，说明窗口特征构造可能有问题，或者 PM2.5 序列本身随机性太强、单靠历史均值无法捕捉突变。这时候不要急着调参，先把时间跨度拉长、数据平滑后再试。此外，验证集上的 RMSE 和训练集上的 RMSE 差异如果很大（比如训练集 R²=0.95、验证集 R²=0.4），就要考虑过拟合——特征维度不多的情况下，过拟合大概率是因为样本量太少。

5. 常见问题与避坑：跑通源码前先看完这五条

这份源码我按「复现 → 调参 → 改数据」的顺序过了一遍，遇到的坑主要集中在数据对齐、学习率稳定性、CSV 编码和模型保存加载几个层面。下面按现象→原因→解决的方式记录五条最高频的问题。

5.1 矩阵乘法报错：shape mismatch

现象：运行 gradient_descent 或预测代码时，numpy 报ValueError: shapes (305,12) and (13,1) not aligned。

原因：训练特征矩阵列数与参数向量行数不一致。通常是拼接 bias 列时只在训练集加了全 1 列，测试集没加；或者滑窗构造特征时窗口长度 window 前后不一致。比如训练时 window=12，测试时却用了 window=11。

解决：在训练和预测两个阶段各打印一次X.shape和theta.shape，确认训练的列数等于 theta 的行数，测试的列数等于 theta 的行数。常见做法是写一个add_bias_column(X)函数，训练和预测都走同一道预处理，避免人工遗漏。

5.2 loss 变成 NaN 或突然爆炸

现象：梯度下降打印的 loss 起初正常，某一次迭代后变成inf或nan，随后一直保持 NaN。

原因：学习率过大导致参数更新步长越过最优区域，loss 在数学上发散，最终溢出为无穷大。另一种可能是数据里存在缺失值，np.nan 参与计算后梯度携带 NaN。

解决：先把 alpha 降到 0.001 试跑，如果还炸，检查 DataFrame 里是否有空值——用pd.isnull(train).sum()逐列检查。PM2.5 序列通常用前一月均值填充缺失值，也可以用插值。还有一个稳健做法是对特征做 min-max 归一化，把数据压缩到 0~1 区间，这样梯度下降的数值稳定性会明显提升；预测完成后再反归一化还原到原始量纲。

5.3 训练集 R² 很高，验证集一塌糊涂

现象：evalu.py 打印训练集 R²=0.92，但用真实未来数据验证只有 0.3 甚至负数。

原因：样本量太少且没有划分验证集，模型「背下」了训练数据的时间波动，对未来数据没有泛化能力。时间序列数据还有一个天然问题是相邻样本高度相关，随机切分训练/测试集会造成信息泄漏，让评估结果虚高。

解决：对时间序列做按时间顺序切分，比如把前 80% 时间点做训练、后 20% 做验证，而不是随机抽样。如果数据只有 12 个月均值，根本撑不起验证集，那就要回到原始按天数据滑窗，把样本量做大。我在评估时一般先打印训练集和验证集的 loss 差值，差值超过 2 倍就该考虑加正则化或减特征。

5.4 model.npy 加载后预测全部是同一个值

现象：predict.csv 里所有预测值几乎相同，变化范围小于 0.01。

原因：特征列的量纲差异过大。PM2.5 若和日期序号（比如 1、2、3…365）拼在同一矩阵，日期序号范围远大于浓度值，梯度下降会优先拟合量纲大的特征，导致其他特征的权重学不到东西。更极端的案例是 theta 初始化全为 0，而学习率太小，1000 次迭代后参数几乎没有移动，输出全部接近 0。

解决：把所有特征全部做标准化（z-score 或 min-max），再进行训练。做完之后哪怕特征里混入日期序号，各列数值范围一致，梯度更新的均衡性也好很多。预测时注意保存标准化器的均值和方差，测试集要使用训练集的统计量做变换，而不是各自重新计算。

5.5 CSV 用 Excel 打开后数据格式错乱

现象：pandas 读入的数组出现字符串类型，np.dot报UFuncTypeError，或数据行数比预期少了一截。

原因：Excel 保存 CSV 时会把长数字改成科学计数法，PM2.5 浓度值如果带一位小数，可能被自动截断或加前缀；中文列名在无 BOM 的 CSV 里被 pandas 默认解析成乱码，导致读错列。

解决：保存数据统一用pd.to_csv(index=False, encoding='utf-8-sig')，utf-8-sig 是中文场景最稳的选择，Excel 能识别、pandas 也能正常读。读取时用dtype={'pm25': float}强制指定列类型，并在读入后打印dtypes做一次检查。如果已经读坏了，用np.loadtxt加delimiter=','按 numpy 方式重新读一遍，通常能绕过 Excel 的格式污染。

6. 把这份源码改成自己的大作业：验证方法、扩展方向和一点建议

拿到一份能跑的源码只是起点。课程大作业最忌讳直接照搬——换个城市、换个预测目标、把模型换成多元线性回归，你的工作量马上就体现出来了。这一章讲怎么在现有代码基础上做稳妥改造，同时给出验证模型可信度的具体方法。

6.1 换成自己的数据：三步替换法

第一步，准备目标城市的空气质量历史数据，至少要包含一列按时间排序的 PM2.5（或 AQI）。第二步，把原始数据塞进 create_dataset 函数生成 arrayx 和 arrayy，注意窗口长度 window 要根据数据粒度调整——月度数据一年 12 行用 window=6 或 12；日度数据 365 行可以用 window=30 或 60，把时间跨度缩短到一个月或两个月。第三步，保持 arrayx.csv、arrayy.csv、x_t.csv 三个文件名不变，直接覆盖原文件，你甚至不需要改主程序就能跑通新数据。

from sklearn.preprocessing import MinMaxScaler scaler = MinMaxScaler() train_scaled = scaler.fit_transform(train[['pm25']]) # 用 train_scaled 构建特征矩阵 X, y = create_dataset(train_scaled.flatten(), window=30) # 做预测后，用 scaler.inverse_transform 还原回真实 PM2.5 值

注意我没有直接用原 12 个月窗口——换成日度数据后，12 天太短，捕捉不到月度季节波动；30 天则对应自然月周期，预测下一天的语义也更合理。这一步的调参空间就在 window 的选择上，建议你分别用 7、30、60 跑一遍，看验证集 RMSE 哪个最小，再写进报告。

6.2 扩展方向：多元线性回归与正则化

如果你想拿高分，把单变量（只用历史 PM2.5）扩展成多元线性回归是最自然的路径。新增特征可以是湿度、风速、气温、气压等气象数据。有了这些列，create_dataset 里就不再是data[i:i+window]一维切片，而是多列同时取窗口、再横向拼接：

def create_multivariate_dataset(df, feature_cols, target_col, window): X, y = [], [] for i in range(len(df) - window): X.append(df[feature_cols].iloc[i:i+window].values.flatten()) y.append(df[target_col].iloc[i+window]) return np.array(X), np.array(y)

flatten() 会把多个特征的窗口值拉平成一维向量：比如 window=7、特征 3 个，每个样本的特征维度就是 21。这种做法能显著提升 R²，但代价是特征矩阵变大，训练样本量要求随之提高。如果你的数据只有几十行，加多元特征很容易过拟合，此时建议引入 L2 正则化（岭回归），在损失函数里加lambda * np.sum(theta[1:] ** 2)，lambda 从 0.01、0.1、1 各试一遍，记录验证集 R² 的变化。

6.3 模型可信度的三条验证习惯

第一个习惯是看 loss 曲线收敛轨迹，不是只看最终 loss 值。把每 100 次迭代的 loss 打印出来画成曲线，如果曲线是平滑下降的，说明学习率合适；如果先降后升，说明 alpha 过大；如果一直锯齿状跳动，需要降低 alpha 或做特征归一化。第二个习惯是算训练误差和验证误差的比值，训练 R² 比验证高 0.3 以上就是一个危险信号，优先减少特征维度而不是增加数据量。第三个习惯是把预测结果还原成真实时间轴画趋势图，看峰值月份是否对应——合肥地区一般冬季 PM2.5 偏高、夏季偏低，如果你的预测曲线趋势方向和真实规律相反，即使数字指标好看也要警觉，多半是数据预处理方向出了问题。