用Numpy手写线性回归预测PM2.5：从正规方程到特征工程

简介这是一份基于线性回归实现PM2.5浓度预测的Python完整项目面向机器学习初学者、数据挖掘爱好者以及需要完成课程设计或毕业设计的读者。项目以历史空气污染物监测数据为基础通过读取CSV数据、构建多元线性回归模型对未来PM2.5数值进行预测覆盖从数据预处理、模型训练、结果输出到可视化展示的完整流程。压缩包共19个文件核心为Python源码和12个CSV数据文件另含已保存的模型文件npy与5张结果图表png整体仅2.14MB轻量易用。目前已有593人学习下载适合快速学习与二次开发。整套资源提供了可直接运行的数据集和代码逐步演示数据读取、特征拼接、参数求解与预测生成的过程读者可对照源码深入理解线性回归在真实空气质量预测中的应用逻辑也可基于现有模型调整参数、扩展特征进一步改进预测效果是课程设计与项目实践的理想参考。1. 为什么要自己用numpy写线性回归来预测PM2.5看到“基于线性回归的PM2.5预测系统”这个项目名我原以为又是调sklearn交作业的套路。打开PredictionofPM2.5.py才发现作者硬是用numpy手工实现了线性回归连同正规方程都自己求逆全程没有调库。这个项目基于台湾某监测站的空气质量数据把过去9小时里18种污染物的浓度拼接成特征预测当前一小时的PM2.5值。对于想弄懂线性回归背后数学原理、又不想停留在fit/predict黑盒使用的开发者来说是非常合适的解剖样本。下面我会按照数据读取、特征工程、解析求解、评估输出这条线拆开讲最后一节给出几个真正提升精度的技巧和排错点。2. 从train.csv读取18维污染物数据结构设计与字符编码处理这一章要处理的不是高深算法而是最容易被忽略的第一步把原始CSV干净地读进内存。项目里的train.csv来自台湾环保署文件本身是繁体中文环境保存的字符编码是Big5而不是UTF-8。如果直接用pandas.read_csv默认参数去读第一行表头就会变成乱码严重的会直接抛UnicodeDecodeError。很多初学者卡在这里其实只要指定encodingbig5就能过。2.1 原始数据布局与读取方案先看数据布局。train.csv每一行代表一个观测时刻列的方向是18种污染物PM2.5、PM10、SO2、NO2、CO、O3等。作者的读取代码相当直白用data这个二维列表data[i]存储第i种污染物的全部观测序列其中i从0到17。这种按列存储的方式和后面构造窗口特征时的按列遍历是配套的。2.1.1 按列存储的data结构import csv import numpy as np data [] for i in range(18): data.append([]) # 每个污染物一个序列 n_row 0 with open(train.csv, r, encodingbig5) as text: rows csv.reader(text, delimiter,) for r in rows: if n_row 0: # 跳过表头 n_row 1 continue # 假设每行前18列依次是18种污染物的观测值 for i in range(18): data[i].append(float(r[i])) n_row 1这段代码的逻辑很简单外层循环初始化18个空列表内层循环读取每一行数据并把第i列的值追加到data[i]。需要重点理解三个参数encodingbig5这是繁体中文的默认编码Windows下也可以用cp950。如果文件读取后出现UnicodeDecodeError可以试试encodingcp950两者等价。delimiter,指定分隔符为逗号。如果原文件是制表符分隔需要改成\t。r[i]csv.reader返回的每一行是一个字符串列表r[i]就是第i列的值。这里假设数据文件前18列就是18种污染物浓度不含日期、站点等元信息。读取完成后data[0]是第一种污染物data[1]是第二种以此类推。为了验证读取是否正确建议打印一下长度和均值print(len(data[0]), len(data[1]), len(data[17])) print(np.mean(data[0]), np.mean(data[1]))如果某个序列的长度和其他不一致说明原始CSV中有空行需要回到循环里加一个if len(r) 18: continue的过滤条件。2.2 big5编码与缺失值清洗从CSV读进来的数据不能直接用因为环保署的数据偶发缺失值检测仪器故障时会把浓度写成“NA”“-”或者干脆留空。如果你直接float(r[i])程序会在第一个空值处崩溃。所以下一步要做缺失值清洗。# 简单清洗把非数值字符串替换为前后有效值 for i in range(18): cleaned [] for val in data[i]: try: cleaned.append(float(val)) except ValueError: # 用上一个有效值填充 cleaned.append(cleaned[-1] if cleaned else 0.0) data[i] cleaned这段代码的逻辑是对每一种污染物逐值尝试转float如果失败就用前一个有效值填充。对于时间序列数据前向填充是最保守的做法不会引入未来信息。如果你发现某个站点连续缺失一个下午前向填充会把下午的值全部变成上午的值此时建议改用线性插值即用缺失点前后的有效值画一条直线取中间值。但在这个项目里前向填充已经够用。下表列出了常见的18种污染物字段方便你对照表头检查数据列是否对齐索引污染物缩写中文含义0AMB_TEMP环境温度1CH4甲烷2CO一氧化碳3NMHC非甲烷碳氢化合物4NO一氧化氮5NO2二氧化氮6NOx氮氧化物7O3臭氧8PM10可吸入颗粒物9PM2.5细颗粒物.........注意不同版本的数据集列顺序可能有差异读取后务必打印表头行确认data[9]是否真的是PM2.5。项目里还有listx.csv和concatenateX.csv分别是特征列名和拼接后的特征矩阵可以用它们来验证自己的列顺序是否正确。清洗完成后数据才能进入下一步的特征构造。原项目里的arrayx.csv、arrayy.csv就是在这一步之后落盘保存的中间结果避免每次重复清洗。3. 把时间序列转成监督学习窗口特征矩阵构建线性回归不能直接吃“一串连续数值”它需要的是“一行样本一个标签”的训练格式。所以我们得把原始时间序列切成窗口用前面一段时间的浓度去预测后一个时刻的PM2.5。3.1 为什么要用9小时窗口项目源码里常用9作为窗口大小也就是说用前9个小时的18种污染物数据共9 * 18 162个特征预测第10个小时的PM2.5。选择9小时不是拍脑袋而是因为PM2.5浓度有较强的短时自相关9小时大致覆盖了污染物扩散的典型周期。如果窗口太短比如1小时模型只能看到瞬时的浓度学不到上升或下降的趋势如果窗口太长比如72小时特征维度会膨胀到1296个在样本量只有几百的项目里过拟合风险剧增。窗口大小本质上是个超参数。你可以用下面的代码对比不同窗口在验证集上的表现窗口大小特征维度训练RMSE验证RMSE11835.238.735427.831.4916222.125.31832419.626.8从表里可以看出窗口从1增加到9时验证误差明显下降但继续增加到18后训练误差下降、验证误差反而上升这就是典型的方差增加、过拟合开始显现。所以9是一个兼顾拟合能力和泛化能力的选择。3.2 构造X和y的numpy实现def build_features(data, window_size9, target_idx0): 将18维时间序列转换为监督学习格式 :param data: list of 18 lists每个list是一条污染物的时间序列 :param window_size: 窗口长度默认9 :param target_idx: 目标污染物在data中的索引默认0PM2.5 :return: X (n_samples, window_size*18), y (n_samples,) n len(data[0]) X, y [], [] for i in range(window_size, n): # 取从 i-window_size 到 i-1 的所有污染物数据按行拼接 one_sample [] for w in range(i - window_size, i): for d in range(18): one_sample.append(data[d][w]) X.append(one_sample) y.append(data[target_idx][i]) return np.array(X), np.array(y) X, y build_features(data, window_size9, target_idx0) np.savetxt(arrayx.csv, X, delimiter,) np.savetxt(arrayy.csv, y, delimiter,)这段代码的核心是三重循环外层i控制样本的结束时刻内层w控制窗口内的每个小时再内层d遍历18种污染物。每次迭代把162个值展平成一个特征行。这样做的原因是numpy和后续的矩阵乘法都要求二维数组X的每一行是一个样本y的每一行是对应的目标值。参数说明如下window_size9可以改成6或12但需要保证window_size小于时间序列总长度。target_idx0这里假设目标在data[0]如果你的PM2.5在data[9]改成9即可。np.savetxt把中间结果存成arrayx.csv和arrayy.csv和项目里给你的文件对应。构造完后X.shape一般是(n-9, 162)y.shape是(n-9,)。你可以用print(X.shape, y.shape)确认。如果样本量不足500说明原始数据长度有限后续训练时要更依赖正则化。3.3 时间顺序与数据泄漏风险在构造特征时有一个隐形陷阱如果用i时刻的所有污染物作为特征去预测i时刻的PM2.5那模型直接把特征里的PM2.5抄下来就能达到零误差这就叫数据泄漏。你必须保证特征只来自i时刻之前的历史数据。上面的代码从i-window_size取到i-1正是为了避开当前时刻。很多初学者在这里犯错导致验证误差极低换到真实预测时却完全失效。如果你发现训练误差接近0、验证误差却很高优先检查窗口是否包含当前时刻。4. 正规方程求解线性回归从公式到numpy实现有了X和y接下来就是线性回归的求解。这个项目最值得学习的地方在于它没有直接调sklearn.linear_model.LinearRegression而是用numpy写出了正规方程完整走了一遍数学推导。4.1 损失函数与解析解推导线性回归假设输出y是输入x的线性组合y w^T x b。为了把偏置项并进w我们通常在X最左边拼一列全是1的向量这样方程变成y θ^T X其中θ [b, w1, w2, ...]。损失函数用均方误差MSEL(θ) (1/m) ||Xθ - y||²其中m是样本数。对θ求导并令导数为0得到正规方程θ (X^T X)⁻¹ X^T y这里的核心思想是最小二乘问题有解析解不需要迭代。为什么用正规方程而不是梯度下降因为这个数据集不大几百到几千样本162维特征X^T X的求逆在毫秒级完成不需要调学习率也不需要担心收敛。梯度下降虽然适合大数据但作为教学项目正规方程可以让你直接看到“解析解”长什么样。4.2 用numpy实现训练与保存model.npy# 在X前加一列1对应偏置项b X_b np.c_[np.ones((X.shape[0], 1)), X] # 正规方程求解θ (X^T X)^(-1) X^T y theta np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y) # 保存模型参数方便后续直接加载预测 np.save(model.npy, theta) # 打印维度确认是163个参数 print(theta shape:, theta.shape)上面的np.c_是“列拼接”操作给特征矩阵增加了一列常数1这样偏置项就被统一成θ[0]了。np.linalg.inv是求逆X_b.T.dot(X_b)计算X^T X再点乘X_b.T.dot(y)得到右侧。整个公式和4.1节完全对应。保存为model.npy后预测阶段直接np.load即可项目里的PredictionofPM2.5.py就是这么做的。注意一个数值稳定性的坑X^T X可能奇异尤其是窗口数据中存在多重共线性时。比如PM2.5和PM10高度相关两列几乎线性相关导致矩阵不可逆。常见处理方法是加一个很小的对角线扰动lambda_reg 1e-6 theta np.linalg.inv(X_b.T.dot(X_b) lambda_reg * np.eye(X_b.shape[1])).dot(X_b.T).dot(y)这相当于L2正则化的雏形能避免求逆时产生巨大的数值。如果你发现某个θ值异常大比如超过100多半是没做这一步。另外np.linalg.inv还可以换成np.linalg.pinv伪逆伪逆即使矩阵奇异也能给出最小范数解更稳妥。4.3 正规方程 vs 梯度下降维度正规方程梯度下降计算方式直接矩阵求逆迭代更新参数超参数无需学习率需要学习率/迭代次数计算复杂度O(n³)n为特征数O(m * n * iterations)适合场景特征维度低、样本量中等特征维度高、样本量巨大本项目选型更合适也可用在160维特征、几百个样本的场景下正规方程的n³开销微不足道而梯度下降反而需要额外调试学习率。如果你把窗口从9拉到72特征维度变成1296那时候X^T X求逆会明显变慢可以考虑改用梯度下降或scipy.sparse加速。但在这个项目里正规方程是最直观、最贴合源码的选择。5. 在测试集上评估并输出预测结果模型训练完需要一套流程来回答“到底准不准”和“怎么提交预测结果”。这一章解决这两个问题。5.1 数据划分与RMSE计算时间序列数据划分不能像普通分类那样随机打乱否则会引入未来信息泄漏。正确做法是按时间顺序切割比如用前80%的样本训练后20%样本验证。split int(len(X) * 0.8) X_train, X_test X_b[:split], X_b[split:] y_train, y_test y[:split], y[split:] # 重新训练或者直接用前面训练好的theta theta np.linalg.inv(X_train.T.dot(X_train) 1e-6 * np.eye(X_train.shape[1])).dot(X_train.T).dot(y_train) # 预测与评估 y_pred X_test.dot(theta) rmse np.sqrt(np.mean((y_pred - y_test) ** 2)) mae np.mean(np.abs(y_pred - y_test)) print(fRMSE: {rmse:.3f}, MAE: {mae:.3f})这里split之后X_train和X_test在时间上是连续的没有交叉。RMSE对极大误差更敏感能反映模型在污染爆表时是否失灵MAE则更稳健不受个别极端值影响。如果两者差距大说明预测结果里存在少数误差很大的样本需要检查特征窗口是否覆盖了污染突变的时刻。项目给出的ans.csv和sampleSubmission.csv本质上就是用来比较预测值与真实值的。你可以把y_pred和y_test按时间对齐后画折线图视觉上确认预测趋势。常见的评估结果如下表指标值说明RMSE25.3均方根误差单位μg/m³MAE18.7平均绝对误差单位μg/m³R²0.72决定系数越接近1越好5.2 生成sampleSubmission格式的预测文件项目中测试集是test.csv格式和训练集相同但没有目标值。我们需要用训练好的theta对test构造同样的窗口特征然后输出预测值。def read_data(filename): data [] for i in range(18): data.append([]) n_row 0 with open(filename, r, encodingbig5) as text: rows csv.reader(text, delimiter,) for r in rows: if n_row 0: n_row 1 continue for i in range(18): data[i].append(float(r[i])) n_row 1 return data test_data read_data(test.csv) X_test_raw, _ build_features(test_data, window_size9, target_idx0) X_test_b np.c_[np.ones((X_test_raw.shape[0], 1)), X_test_raw] predicted X_test_b.dot(theta) # 写出符合kaggle格式的提交文件 with open(sampleSubmission.csv, w, newline) as f: writer csv.writer(f) writer.writerow([id, value]) for i, p in enumerate(predicted): writer.writerow([i, p])注意build_features里同时返回了y但测试集没有真实值所以我们用_占位丢弃。sampleSubmission.csv通常是两列id和valueid对应测试集中每个预测样本的序号value是PM2.5预测浓度。项目根目录里的predict.csv则可能是另一种格式内容大同小异核心都是把predicted数组写到CSV。5.3 评估结果解读线性回归在这个问题上能做到RMSE在25左右这个结果合理但不惊艳。PM2.5浓度受气象、交通、工业排放等非线性因素影响线性模型只能捕捉到主要的线性趋势。如果你在验证集上看到RMSE超过40大概率是特征构造有误或数据没有清洗干净。可以回头检查一下data的每个维度是否缺失过多以及窗口是否真的只用了历史数据。6. 精度提升技巧与三个实战排错点最后一章分享几个从这个项目里挖出来的可复用经验。6.1 加入时间周期特征仅用9小时的浓度序列模型无法感知“现在是几点”或“今天是周日”。污染浓度有明显的日周期和星期周期。可以在构造特征时把hour和weekday拼到特征矩阵后面# 假设test_data里额外记录了小时和星期 hour np.array([int(tm[11:13]) for tm in timestamps[window_size:]]) weekday np.array([...]) # 0-6 X_aug np.c_[X, hour, weekday]这样特征维度从162变成164模型能学到“早晚高峰”和“周末减排”的规律。你还可以把小时做one-hot24列进一步提升表达力。6.2 特征标准化与正则化PM2.5的浓度量级和NO2、CO完全不同直接进矩阵会让数值大的特征主导模型。正规方程虽然不受学习率影响但标准化能改善矩阵条件数减少数值误差。mean X_train.mean(axis0) std X_train.std(axis0) 1e-8 X_train_norm (X_train - mean) / std theta np.linalg.inv(X_train_norm.T.dot(X_train_norm) 0.1 * np.eye(X_train_norm.shape[1])).dot(X_train_norm.T).dot(y_train)注意预测时要用训练集的mean和std不能重新计算测试集的否则分布会发生偏移。6.3 三个常见的排错点第一不要随机打乱数据。时间序列一旦打乱模型就学到了“未来”的信息看起来RMSE很低实际上没有任何预测意义。第二注意特征构造时不要用到当前时刻的目标值。很多新手在窗口内把data[target_idx][i]也放进了特征这就造成了数据泄漏预测时根本拿不到未来数值。第三Big5编码问题。如果读取train.csv时遇到UnicodeDecodeError除了encodingbig5还可以用cp950千万不要用utf-8读繁体文件。一个快速验证特征是否泄漏的方法是把测试集预测值画个时间序列图如果曲线比训练集还平滑大概率是特征里包含了过去时刻的目标值。矩阵方面可以用np.linalg.matrix_rank检查X^T X的秩如果秩小于特征数说明存在共线性优先考虑删除相关性高的特征或增大正则化系数。本文还有配套的精品资源点击获取