☰
基于线性回归的PM2.5预测实战:从数据预处理到模型评估
2026/10/9 14:39:00 网站建设 项目流程

简介:这是一份面向机器学习初学者的Python课程大作业源码包,选择合肥地区过去一年的PM2.5月度数据作为样本,完整实现基于线性回归的空气质量预测。项目不仅包含数据读取与清洗、梯度下降公式推导与代码实现、矩阵模型构建,还提供了训练好的模型文件和预测结果示例,可直接运行或对照学习。包内共21个文件,主要有3个Python脚本、12个CSV数据集、1个npy模型文件、若干图片和说明文档,整体约2.58MB,目录中训练集、测试集、中间数组及结果文件划分明确,方便按步骤排查与调参。目前已有280人学习/下载。通过源码可以掌握线性回归在真实数据集上的应用流程,学会利用梯度下降更新参数、用矩阵运算加速计算,并结合可视化图片、结果示例与项目文档理解课程设计中的关键逻辑,适合作为课程作业参考或算法入门实践。

1. 机器学习大作业-基于线性回归的PM2.5预测源码:能复现到什么程度

如果你正在为一门机器学习课程设计找现成项目,又不想从零调参,"基于线性回归的PM2.5预测"这种题目看起来朴素,实际坑相当多。这个源码包正是拿合肥地区过去一段时间(比如过去一年每月的平均值)的 PM2.5 数据,训练一个线性回归模型,预测今年某个月的空气质量值。它把训练、预测、评估、画图四件事都做完了,还留下 train.csv、test.csv、模型权重、中间特征矩阵等一整套文件,不是那种只有一个主文件的半成品。适合正在做课程大作业、需要复现并在答辩时讲清楚原理的人——先跑通一遍,再换成自己所在城市的数据就能交差。

2. 数据组织与特征拼接:先把 train.csv 和 test.csv 的格式对齐

2.1 先认识这几个 csv:哪些是输入,哪些是中间产物

解压后是一个 Machine-Learning-master 目录,README.md 和课程样例给了基本使用说明,但真正干活的是下面这批文件。我建议拿到手第一件事不是直接跑 PredictionofPM2.5.py,而是把 csv 文件按角色分清楚,避免后面读错文件。

文件角色作用
train.csv原始训练输入合肥地区历史 PM2.5 数据,可能包含日期和多个污染物列
test.csv / testdata.csv测试输入与 train.csv 同格式,用于构建测试特征
sampleSubmission.csv提交样例输出格式参照,最后 predict.csv 可以仿照它写
arrayx.csv特征工程中间产物滑动窗口切出来的特征矩阵 X
arrayy.csv特征工程中间产物每个窗口对应的真实值 y
x_t.csv / concatenateX.csv特征矩阵调试产物转置或拼接了偏置列的最终矩阵
predict.csv预测结果模型在测试集上的输出
ans.csv评估对照测试集真实值,evalu.py 靠它算误差
s_gra.csv绘图数据真实值、预测值、月份的对照表
model.npy模型权重训练完的 theta,预测时直接加载
image.png / demo.jpg / a th.jpg可视化结果损失曲线或预测对比图

这些文件里真正要亲手改的是 train.csv 和 test.csv,其余中间产物都是脚本生成的。第一批踩坑往往发生在读数据阶段:train.csv 的列名和 README 里写的不完全一致,或者日期列解析失败。

先做一步读取确认:

import pandas as pd raw = pd.read_csv('train.csv') print(raw.columns.tolist()) print(raw.head())

逻辑说明:先看列名,再决定后续按哪一列解析日期、哪一列作为 PM2.5 数值。很多课程数据集的列名带空格或大小写混用,直接 groupby 会报 KeyError。参数说明:read_csv 的 encoding 参数在遇到中文列名时常用 utf-8-sig,否则 Windows 下容易乱码。

合肥的数据如果是小时级监测记录,需要先聚合到月。摘要里写的是"过去一年每个月的平均值",所以这个聚合步骤是核心前置操作。

def to_monthly(df, value_col='pm2.5'): df = df.copy() df['date'] = pd.to_datetime(df['date']) df['month'] = df['date'].dt.to_period('M') return df.groupby('month')[value_col].mean().reset_index() monthly = to_monthly(raw) print(monthly)

逻辑说明:to_datetime 负责把字符串日期转成时间类型,dt.to_period('M') 把所有日期归到所在月份,最后 groupby 按月求平均,得到类似"2024-01、2024-02..."的月均值序列。参数说明:value_col 指定目标列,如果原始列叫 PM2.5 而不是 pm2.5,这里要手动对齐;groupby 之后会丢掉其他污染物信息,这个作业只需要 PM2.5 单变量回归,所以没问题。

2.2 特征矩阵怎么拼:从 arrayx.csv 到 concatenateX.csv

线性回归做时间序列预测,核心思路是用过去连续 N 个月的月均值作为特征,预测下一个月。这份资源里 window 一般取 6,也就是拿前半年拟合后半年。

import numpy as np window = 6 X, y = [], [] for i in range(window, len(monthly)): X.append(monthly['pm2.5'].values[i - window:i]) y.append(monthly['pm2.5'].values[i]) X = np.array(X).reshape(len(X), -1) y = np.array(y).reshape(-1, 1) np.savetxt('arrayx.csv', X, delimiter=',') np.savetxt('arrayy.csv', y, delimiter=',')

逻辑说明:循环从第 window 行开始,每次取前 6 个月的 PM2.5 值作为一条样本,第 7 个月的值作为标签。数组最终形状是 (样本数, 6),每一行就是一个滑动窗口。参数说明:window 太小模型学不到趋势,太大合肥这种一年只有 12 个月均值的数据样本量会急剧减少,6 是一个比较平衡的值。样本量不足时,也可以把 window 降到 3,代价是特征信息变少。

这里的 listx.csv 我一般理解为窗口的起始月份编号,用来在答辩时说明"第 i 行特征对应的是哪 6 个月"。它不影响训练,但排查特征错位时很有用。x_t.csv 则是一次调试时把特征矩阵转置后保存的视图,方便核对矩阵维度,主流程不依赖它。

真正交给模型的是 concatenateX.csv,它比 arrayx.csv 多了一列全 1:

Xb = np.hstack([np.ones((X.shape[0], 1)), X]) np.savetxt('concatenateX.csv', Xb, delimiter=',')

逻辑说明:线性回归的矩阵形式 h_theta(X) = X·theta 里,theta 的第一项对应截距 bias。如果不拼这列全 1,模型只能拟合过原点的直线,PM2.5 月均值永远不会是 0,预测必然系统性偏低。参数说明:np.hstack 是横着拼列,要求两个数组行数一致,这也是为什么不直接改 arrayx.csv,而是另存一个 concatenateX.csv 的原因——保留中间产物便于反复检查。

3. 核心实现:PredictionofPM2.5.py 里的矩阵模型与梯度下降

3.1 模型形式与两个求解方式

这份作业的模型是标准多元线性回归:h_theta(X) = X·theta。X 是上一步拼好的特征矩阵,theta 是待学习的权重向量,输出是 PM2.5 月均值。梯度下降公式在 README 和课程样例里都给了,核心是矩阵形式批量更新:

theta := theta - (alpha / m) * X.T · (X·theta - y)

这个公式里 X.T·(X·theta - y) 是一次性算出所有样本的梯度方向,再取平均。为什么要用矩阵形式而不是写 for 循环逐样本更新?因为矩阵运算一次把所有样本的误差聚合了,训练 2000 轮也就 2000 次矩阵乘法,而逐样本更新在月均值这种样本数只有几十的数据集上虽然也能收敛,但代码啰嗦且容易在答辩时被追问"为什么不用向量化"。

正规方程是另一个闭式解,theta = (X.T·X)^(-1)·X.T·y,一步出结果。但作业明确要求梯度下降公式,所以主脚本走迭代路线,正规方程可以作为对照实验。两者结果应当非常接近,如果差异过大,说明特征没对齐或标准化处理不一致。

梯度下降的超参数一般这样设:

参数含义常见取值说明
alpha学习率0.01PM2.5 数值在几十到几百,标准化后 0.01 起步安全
iters迭代次数2000 到 5000看损失曲线是否已收敛
window特征窗口6过去 6 个月预测下一个月
m样本数由数据决定12 个月年均值时样本数只有 6 个左右

3.2 训练与保存主流程:标准化、迭代、落盘 model.npy

PredictionofPM2.5.py 的主流程可以拆成四步:读 concatenateX.csv 和 arrayy.csv、StandardScaler 标准化、梯度下降训练、np.save 保存 theta。标准化这一步非常关键,合肥 PM2.5 月均值秋冬能到 120+,夏季可能只有 30,特征列之间数值差异大,梯度下降在这样的尺度下很容易震荡。

import numpy as np import pandas as pd from sklearn.preprocessing import StandardScaler X = pd.read_csv('concatenateX.csv', header=None).values y = pd.read_csv('arrayy.csv', header=None).values scaler = StandardScaler() X_scaled = scaler.fit_transform(X[:, 1:]) X_scaled = np.hstack([np.ones((X_scaled.shape[0], 1)), X_scaled]) theta = np.zeros((X_scaled.shape[1], 1)) alpha, iters = 0.01, 3000 def gradient_descent(X, y, theta, alpha, iters): m = len(y) cost_history = [] for i in range(iters): h = X.dot(theta) errors = h - y theta = theta - (alpha / m) * X.T.dot(errors) cost = float(np.mean(errors ** 2) / 2) cost_history.append(cost) if len(cost_history) > 1 and abs(cost_history[-2] - cost) < 1e-8: break return theta, cost_history theta_final, costs = gradient_descent(X_scaled, y, theta, alpha, iters) np.save('model.npy', theta_final)

逻辑说明:fit_transform 只对真实特征列做标准化,bias 那列全是 1,不需要也不应该被标准化。梯度下降内层循环先算预测值 h,再算残差 errors,然后用 X.T.dot(errors) 得到梯度,最后更新 theta。cost_history 记录每轮损失,用于判断收敛。参数说明:m 是样本总数,学习率除以 m 是为了把梯度平均到每个样本;1e-8 的阈值是提前停止条件,损失变化小于这个值就认为收敛,可以手动调小但没必要,3000 轮在这个数据量下通常一两秒就跑完。

预测脚本加载 model.npy 时有一个容易忽略的点:model.npy 里只存了 theta,没有存 scaler 的 mean 和 std。所以预测时要么重新读训练集做一遍 fit,要么把均值方差也存下来。这份资源里预测脚本的做法是重新加载 train.csv 重建特征再做同样的标准化,实际跑的时候我建议直接改成保存三个数组到同一个 npz 里,省得每次重新 fit。

注意:np.save 保存的是一维数组时,加载后要确认 shape 是否为 (特征数, 1),否则后续矩阵乘法会广播出形状错误。

4. 评估与可视化:evalu.py 和 s_gra.csv 怎么配合用

4.1 评估指标按什么口径算

训练完不能只看训练集损失小就交差,测试集才是老师真正看的东西。这份资源里的 evalu.py 读 predict.csv 和 ans.csv,逐行对比预测值和真实值。评估指标建议算三个:RMSE、MAE、R²。

import numpy as np import pandas as pd pred = pd.read_csv('predict.csv')['pm2.5'].values true = pd.read_csv('ans.csv')['pm2.5'].values rmse = np.sqrt(np.mean((pred - true) ** 2)) mae = np.mean(np.abs(pred - true)) r2 = 1 - np.sum((true - pred) ** 2) / np.sum((true - np.mean(true)) ** 2) print(f'RMSE={rmse:.2f} MAE={mae:.2f} R2={r2:.3f}')

逻辑说明:RMSE 对大误差敏感,预测值里如果有一个月和真实值差了 40 µg/m³,RMSE 会被明显拉高;MAE 更平均地反映整体偏差;R² 表示模型解释了真实值方差的多少,越接近 1 越好。参数说明:真实值和预测值里如果存在 NaN,pandas 默认在求和时会变成 NaN,所以评估前最好加 np.nan_to_num 或 dropna。合肥月均值量级一般在 30 到 120 之间,RMSE 在 10 左右已经算不错,R² 在小样本上容易虚高,不要只看 R²。

4.2 把真实值和预测值画到同一张图

s_gra.csv 是绘图直接的数据源,里面通常有 month、true、pred 三列。训练完模型后,用测试集月份生成预测值,再把真实值和预测值拼进同一个 DataFrame,最后存成 s_gra.csv 并画图。

import matplotlib.pyplot as plt import pandas as pd df = pd.read_csv('s_gra.csv') plt.figure(figsize=(10, 4)) plt.plot(df['month'], df['true'], label='true', marker='o') plt.plot(df['month'], df['pred'], label='pred', marker='s') plt.legend() plt.savefig('demo.jpg', dpi=200) plt.close()

逻辑说明:真实值用实线圆点,预测值用方块线,两张图叠在一起能直观看出哪个月预测偏大、哪个月偏小。demo.jpg 应该就是这一步的产物,image.png 和 a th.jpg 可能是不同窗口或不同学习率下的对比图。参数说明:dpi=200 保证图片在论文里插进去不模糊;列名如果和上面代码不一致,先用 df.columns.tolist() 看实际列名,再替换 df['month'] 里的字段名。

s_gra.csv 的另一个用途是存档。答辩时老师问"你哪个月预测最差",直接打开 csv 按残差排序就能回答。我一般会额外加一列 abs_error,方便按误差大小过滤。

df['abs_error'] = (df['true'] - df['pred']).abs() print(df.sort_values('abs_error', ascending=False).head())

逻辑说明:按绝对误差降序排,排在最前面的就是模型表现最差的两个月,通常对应秋冬季节的突变天。参数说明:ascending=False 是降序,head(3) 取前三行,用来在答辩材料里写"模型在污染突变月份误差较大"这种结论。

5. 避坑指南:合肥PM2.5预测里最容易翻车的五个细节

5.1 损失不降反升:指数爆炸和 NaN 都是同一个原因

现象:gradient_descent 循环跑完,print(cost_history) 发现损失不是下降而是每次都变大,甚至跑到几十轮后出现 NaN。原因:特征列没标准化,PM2.5 数值在 120 时 X.T.dot(errors) 的结果量级会非常大,alpha 虽然设了 0.01,但乘上去之后 theta 直接飞掉。解决:先做 StandardScaler,alpha 从 0.001 开始试,并且每 50 轮打印一次损失:

alpha = 0.001 for i in range(iters): ... if i % 50 == 0: print(f'iter {i}, cost {cost:.4f}')

逻辑说明:alpha 太小收敛慢,alpha 太大直接震荡,0.001 到 0.01 之间是标准化后的安全区间。参数说明:打印频率 50 轮一打印足够看到趋势,不用每轮都打。

5.2 时间序列不能随机 split:不小心就偷看了未来

现象:作业里用了 train_test_split(X, y, random_state=42),训练集 RMSE 好得离谱,但换成测试集后预测曲线明显滞后。原因:时间序列数据一旦随机划分,后半段的样本会混进训练集,模型提前学到了"未来"的走势,测试集反而只留下随机噪声。解决:按时间顺序划分,train_test_split 里强制 shuffle=False:

from sklearn.model_selection import train_test_split X_tr, X_te, y_tr, y_te = train_test_split( X, y, test_size=0.3, shuffle=False, random_state=0 )

逻辑说明:shuffle=False 保持原始顺序,前面的月份进训练集,后面的月份进测试集,模拟真实预测场景。参数说明:random_state 固定成 0 是为了每次跑结果一致;如果数据只有 12 个月均值,test_size 不要取太大,留 3 个月做验证差不多。

5.3 model.npy 的维度对不上:忘了 bias 那一列

现象:训练时读的是 concatenateX.csv,特征矩阵有 7 列(1 列 bias + 6 列特征);预测时只读了 6 列特征,结果 np.load('model.npy') 之后 dot 报错,或者没报错但所有预测值都整体偏移一个常数。原因:theta 的第一个分量是截距 bias,预测脚本必须给输入矩阵同样拼一列全 1。解决:把训练时的偏置拼接逻辑复制到预测脚本:

X_test = pd.read_csv('test.csv', header=None).values X_test_scaled = scaler.transform(X_test) X_test_scaled = np.hstack([np.ones((X_test_scaled.shape[0], 1)), X_test_scaled]) pred = X_test_scaled.dot(theta_final)

逻辑说明:预测端的特征构造必须和训练端完全一致,包括窗口大小、特征顺序、是否加 bias 列。参数说明:scaler 要用训练时那个 fit 过的 scaler,不能在测试集上重新 fit,否则分布就变了。

5.4 "0值"到底是 0 还是缺测:噪声数据会拖低整条曲线

现象:模型在冬季月份预测整体偏低,残差图里有一两个点特别尖。原因:PM2.5 监测仪器在故障或断电时会输出 0,这个 0 不是真实浓度,而是机器学习的噪声数据。如果直接用 0 值参与训练,模型会以为那几个月空气质量特别好,把权重往低拉。解决:把等于 0 的值先替换成 NaN,再插值:

monthly['pm2.5'] = monthly['pm2.5'].replace(0, np.nan) monthly['pm2.5'] = monthly['pm2.5'].interpolate()

逻辑说明:replace(0, np.nan) 把噪声标记成缺失,interpolate 用线性插值填补。参数说明:interpolate 默认线性法,对月均值这种低频序列足够;如果连续多个月都是 0,说明那段时间数据质量太差,直接剔除比填充更安全。

5.5 中间 csv 的列名是最好的排查抓手

现象:读 arrayx.csv 时,pandas 自动把列名设成 0 到 5;读 concatenateX.csv 时又把第一行当成了列名,导致特征矩阵少一行或维度对不上。原因:np.savetxt 保存时没有写表头,而 pd.read_csv 默认认为第一行就是列名。解决:读中间 csv 一律显式指定 header=None,再用数值下标取列:

X = pd.read_csv('arrayx.csv', header=None).values y = pd.read_csv('arrayy.csv', header=None).values

逻辑说明:header=None 告诉 pandas 所有行都是数据,列名用 0 开始的整数替代。参数说明:.values 会把 DataFrame 转成 ndarray,后续矩阵运算需要的是 ndarray 而不是 DataFrame,这个转换别省。

6. 进阶:把模型权重拿出来做特征排序和提前收敛

6.1 标准化后的系数绝对值就是特征重要度

线性回归出了名的可解释性强。模型训练完,model.npy 里存着的 theta 直接反映了每个输入特征对 PM2.5 的贡献方向。特征全部标准化后,权重绝对值越大,说明该月份对预测目标的影响越强。

theta = np.load('model.npy').flatten() weights = theta[1:] rank = np.argsort(np.abs(weights))[::-1] for i, idx in enumerate(rank): print(f'前第{i+1}重要: 第{idx+1}个月前值, 权重 {weights[idx]:.3f}')

逻辑说明:theta[1:] 去掉 bias,因为截距只负责整体平移,不参与特征排序。np.argsort 返回排序后的索引,[::-1] 翻转成从大到小。输出结果可以直接写进大作业的"特征分析"章节,解释为什么春秋季节的 PM2.5 更容易受前一个月影响。参数说明:这里比较的前提是训练时做了 StandardScaler,否则不同月份数据的数值尺度不同,权重不可比。

6.2 用损失曲线判断要不要加大迭代次数

gradient_descent 返回的 cost_history 不只是用来确认收敛,还能判断当前数据量下模型是否欠拟合。把曲线画出来:

import matplotlib.pyplot as plt plt.figure(figsize=(8, 3)) plt.plot(costs) plt.xlabel('iterations') plt.ylabel('cost') plt.savefig('image.png', dpi=150) plt.close()

逻辑说明:如果曲线在最后几十轮还在明显下降,说明 iters 不够,加迭代;如果曲线快速收敛后平坦,说明 3000 轮足够。image.png 在原始文件里就是一张这类图,我推测就是某次训练保存下来的损失曲线。参数说明:costs 是梯度下降函数返回的列表,长度等于实际迭代轮数,提前停止会让它小于 iters。

从那以后我每次提交大作业前,都会强制走一遍完整流程:训练模型、跑 evalu.py、看 s_gra.csv 画出的曲线,确认 RMSE 没有突然恶化、预测值和真实值没有整体错位,再导出 predict.csv。这个习惯帮我挡过至少两次因为忘加 bias 列而导致的翻车。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询