线性回归预测PM2.5：从数据预处理到梯度下降的完整复现

最近我把一门机器学习公开课的2020年版作业1翻出来重新做了一遍题目就是预测PM2.5的值要求用线性回归完成。这个作业在圈子里流传很广很多人第一次接触机器学习都是从这道题开始的。它虽然叫“线性回归”但真正做起来会发现数据清洗、特征构造、训练验证的坑一个都不少。这里把我完整复现的思路、代码和踩坑记录整理出来希望能帮你少走点弯路。1. 这个作业到底在考什么线性回归不止是“套公式”1.1 作业任务的一句话版给定某城市空气质量监测站记录的逐小时数据每天有24个小时每小时观测到18个特征其中包括PM2.5、PM10、一氧化氮、二氧化氮、臭氧等污染物浓度以及温度、湿度、风速等气象指标。要求是用连续9个小时的历史观测值预测第10个小时的PM2.5浓度。听起来很简单就是找一条直线或者一个超平面去拟合数据。但作业真正的考点有三层第一层会不会把原始数据整理成模型能吃的训练样本。这是基础门槛。第二层能不能理解线性回归的两种求解方式也就是闭式解和梯度下降并解释它们的区别。第三层能不能基于验证集判断模型效果而不是只把训练集loss压得很低。很多初学者上来就调库sklearn.linear_model.LinearRegression一行代码跑完然后输出答案交作业。这样当然也能过但失去了这个作业最大的价值。我下面会完全用NumPy手写线性回归不依赖机器学习框架这样能看清每个参数在更新的时候到底发生了什么。1.2 这个数据集有哪些隐藏的细节我最初拿到数据的时候第一反应是“这不就是个表格吗”但打开一看就开始头疼了。先看结构原始数据是CSV文件按月份组织每个文件里有18行特征指标每一列是一个小时。对你没有看错常规表格是一行一个样本但这个数据集是“一行一个特征指标一列一个时间点”。18行对应18项监测指标每一行有24个小时的观测值所以一个月大约有20天的训练数据时这个文件的维度就是18x480左右。这种存储格式直接照搬进模型是不行的。做预测时我们希望训练样本长这样每一行是“过去9小时所有特征的观测值”也就是一个长度为9乘以18的向量然后塞进线性回归里。所以第一步必须做转置和窗口切片。另一个容易忽略的细节是数据里有NR这种字符型缺失值标记。我查了一下NR在原始资料里代表该项仪器当时没有读数。如果用pandas直接读取这一列会被读成字符串后续做矩阵乘法直接报错。清洗时必须把它替换成数值我习惯用0填充也可以用前后时间点的均值填充。两种做法我都试过对最终预测准确率影响不大但处理方式必须写在代码里不然第二天自己都看不懂。1.3 适合谁按照这个思路做如果你是第一次接触机器学习我建议你不要跳过数据预处理部分直接看训练代码是没有意义的。这篇文章更适合下面几类人刚学完线性代数想看看矩阵运算在真实数据上怎么落地的人。正在做公开课作业但卡在“从原始CSV到训练样本”这一步的人。已经交过作业但回头想优化预测精度想知道除了调库还能怎么做的人。2. 从原始监测表到训练样本数据预处理是重头戏2.1 把“18行x很多列”变成连续时间序列我先说我的做法用pandas读取每个月的CSV把所有月份纵向拼接成一个大表。注意这里的“纵向拼接”指的是时间维度的拼接不是特征维度的拼接。处理思路分三步。第一步把每个月的文件读进来丢弃不需要的指标列。原始数据里有18个特征但有的指标缺失率特别高比如某个监测项几乎一半以上都是NR。这种特征直接放进模型不仅提供不了有效信息还会增加训练时间。我第一版保留了全部18个特征第二版只保留PM2.5、PM10、NO2、SO2、O3和温度、湿度、风速这8个效果差别不大但训练速度快了一些。第二步把读取到的数据转置。转置之后DataFrame的每一行就变成了某个小时的观测值列就是18个特征。因为CSV是按月组织的转置后还要把多个月份的数据按时间顺序拼接起来。第三步把拼接后的DataFrame显式构造出“日期时间”这个维度。这个很重要。原文件里没有明确的时间索引只是一列一列地排下来如果你不自己维护一个小时的索引后面做滑窗的时候就很容易把不同日期的数据混在一起。2.2 用滑窗构造“9小时预测下一小时”的样本假设我们已经有了一条连续的时间序列长度为N每条记录有18个特征。我们要建立的映射关系是(t-8, t-7, ..., t时刻的特征) (t1时刻的PM2.5)也就是说用当前时刻以及过去8个小时的所有特征预测下一个小时的PM2.5。这里有人会疑惑为什么不是用过去9个小时预测当前时刻如果是给定第1到第9小时的数据预测第10小时那么用于构造样本的规律应该是——特征窗口从第i到第i8小时标签是第i9小时的PM2.5。两种写法本质等价但注意别把对应关系搞反。滑窗具体操作是这样的设定窗口长度window_size 9。从时间序列的第0行开始取第0行到第8行这9行展平后是一个长度为9*18162的向量。标签取第9行的PM2.5值。然后把窗口向右移动一格也就是第1行到第9行作为新的特征向量第10行的PM2.5作为标签。一直滑到序列末尾。滑窗会得到N-9个样本。比如连续500个小时的数据能得到491个样本。这里的边界问题是天然存在的序列中最后几个小时的记录无法构造出“未来一小时”因为没有了未来数据所以在任务给定的测试集里通常已经帮你切好了样本不需要你处理滑窗的首尾。我踩过的一个坑是原始月份之间是断开的。比如一个月最后一天晚上23点的数据和下个月第一天0点的数据从时间上看是连续的但它们属于不同的文件读取批次。如果直接把它们当成连续序列做滑窗就可能构造出跨月边界的、物理上不存在的训练样本。我当时的解决方案是在每个月份文件拼接时记录边界索引只在同一个月的数据内部滑窗跨月边界处断开。2.3 归一化的时机和方式滑窗构造完样本矩阵后特征向量里每个维度的数值范围差异非常大。PM2.5浓度可能只有个位数到一两百总悬浮颗粒物可能从几微克到几百微克温度可能是十几度到三十几度风速也可能是0.1到10以上。如果直接把原始数值丢给梯度下降loss曲线会非常难看可能来回震荡几千轮都收敛不了。归一化的方式我选的是标准分数也叫Z-scorez (x - mean) / std这里的均值和标准差只应该在训练数据上计算。我见过有人偷懒对全部数据一起算这在作业题里不是不行但如果希望模型评估严谨些必须只用训练集拟合归一化参数再用同样参数去归一化测试集。否则相当于让模型提前偷看了测试集的数据分布验证结论会有水分。归一化有个容易忽略的点多个特征共享一个均值还是每个特征单独算正确答案是每个特征单独算。因为我们展平后的向量维度是162也就是9个时间步乘以18个特征每一维都有自己独立的均值和标准差不能用整个样本矩阵的全局均值和全局标准差。处理之后样本矩阵每一列的数值都大致落在[-3, 3]的范围内这才适合接下来做矩阵运算。3. 线性回归的两种打开方式闭式解与梯度下降3.1 目标函数先说清楚线性回归假设y Xw b其中y是预测的PM2.5浓度X是一个行向量包含162个特征w是对应的权重向量b是偏置。为了表达简洁通常把偏置合并进权重向量在X后面拼一个常数1相应地w多一个维度。这个技巧叫“增广特征”。损失函数用均方误差L (1/(2m)) * sum((Xw - y)^2)这里m是样本个数。前面这个1/2不是必须的纯粹是为了求导时把2抵消掉很多教材喜欢这么写。实际运算时有没有这个系数只影响损失数值的绝对大小不影响最优解。3.2 闭式解直接算一步到位当特征数量不是特别大且矩阵可逆时可以直接用正规方程求解w* (X^T X)^(-1) X^T y这个公式的推导过程不复杂对损失函数求导让导数为0整理以后就是上面这个式子。我在实际做的时候遇到的最大问题不是公式记不住而是矩阵太大算不动。如果我用全量训练数据构造出的X维度是(m, 162)那么X^T X的维度就是(162, 162)求逆的代价其实很小。真正的问题是如果样本量m达到几十万甚至上百万构造X^T X的过程需要做一次超大矩阵乘法内存吃紧。我在自己的机器上试了一下训练数据大约30万条样本时X^T X算得还算流畅但如果样本量翻倍时间会非常感人。一个更实际的问题是正则化。如果特征之间存在多重共线性或者特征数接近样本数X^T X可能不可逆。解决办法是加一个很小的对角扰动w* (X^T X λI)^(-1) X^T y这个λ就是L2正则化系数作业里一般不会强制要求但加上会稳很多也不容易过拟合。3.3 梯度下降理解参数怎么动起来闭式解适合特征维度不高、一次算完的场景。但作业明确要求用梯度下降原因是它更能体现机器学习迭代优化的思想。梯度下降的更新公式w - w - η * (1/m) * X^T (Xw - y)这里η是学习率。我看到很多新手把公式抄下来就开跑然后损失在某个值附近疯狂震荡。原因多半是学习率没调。我做了几组对比实验学习率训练loss变化情况结论0.01下降缓慢5000轮还没到收敛太低耗时间0.1稳步下降约1000轮收敛适合该数据集0.5前100轮下降快之后震荡略高需要调低1.0loss直接爆炸太高发散值得注意的是这个学习率是否合适和归一化强相关。如果不归一化直接用0.1的学习率大概率也会发散因为不同特征的梯度量级差异太大。我在实验里先归一化再训练学习率0.1就比较好用。梯度下降过程中还有一个容易被忽视的点要不要保存历史loss。我强烈建议每隔一定轮数打印一次loss画成曲线。如果loss曲线出现“先降后升再降”的诡异形态说明学习率或者数据划分有bug早点发现能省出大量排错时间。3.4 为什么很多人第一版只挑PM2.5作为特征题目虽然给了18个指标但很多公开解法只取PM2.5这一个指标来预测未来一小时PM2.5也就是把162维特征压缩成9维。这是最简化的线性回归示例优点是训练快、代码简单缺点是模型表达能力有限。我在复现时两种都做了只使用PM2.5过去9小时数据特征维度9。使用全部18个指标过去9小时数据特征维度162。对比下来全部特征在验证集上的RMSE大概比只用PM2.5低15%到25%。原因很好理解空气质量不只是PM2.5自身演化的结果还受气象条件和其他污染物影响。比如风速提升有利于颗粒物扩散湿度升高可能导致二次颗粒物生成。这些信息不在PM2.5的历史序列里只能靠其他指标提供。所以我的建议是作业验证阶段可以先跑通9维特征看全流程没问题后再扩展到162维特征这本身就是一次超实用的实验对比。4. 手写NumPy实现从零训练一个线性回归模型4.1 数据读取与样本构造代码这一段代码是把前面第2节的数据处理思路落地的完整实现。我用的是Python加NumPy不需要其他机器学习库。import numpy as np import pandas as pd # 读取某个月的原始CSV def read_month_data(file_path): df pd.read_csv(file_path, encodingutf-8) # 原始文件第一列是时间指标名称后续每一列是一个小时 # 转置让每一行成为一个小时的观测 df df.transpose() # 第一行变成表头 df.columns df.iloc[0] df df.drop(df.index[0]) # 把NR替换成0.0 df df.replace(NR, 0.0) return df.astype(float) # 滑窗构造样本 def create_samples(df, window_size9): all_x [] all_y [] data df.values num_hours, num_features data.shape for i in range(num_hours - window_size): feature_window data[i : i window_size, :] label data[i window_size, 3] # 假设PM2.5是第4列索引为3 all_x.append(feature_window.flatten()) all_y.append(label) return np.array(all_x), np.array(all_y)我这里假设PM2.5在原始特征里索引为3实际使用时要先打印出列名确认一下。不同的CSV版本里列顺序可能有差异这是个必须手动检查的点不能盲目相信列名顺序。4.2 归一化与划分训练验证集构造完X和y之后要把数据集分成训练集和验证集。我选择的划分方式是按时间顺序前80%的数据作为训练集后20%作为验证集。这样做的理由是空气污染数据具有时间相关性随机打乱再划分虽然也能用但会乐观地高估模型效果因为模型见过邻近时刻的数据。def split_and_normalize(X, y, train_ratio0.8): split_idx int(len(X) * train_ratio) X_train, X_test X[:split_idx], X[split_idx:] y_train, y_test y[:split_idx], y[split_idx:] mean X_train.mean(axis0) std X_train.std(axis0) std[std 0] 1.0 # 避免除零 X_train (X_train - mean) / std X_test (X_test - mean) / std # 给特征加一列常数1对应偏置项 X_train np.concatenate([X_train, np.ones((X_train.shape[0], 1))], axis1) X_test np.concatenate([X_test, np.ones((X_test.shape[0], 1))], axis1) return X_train, X_test, y_train, y_test这里对std为0的列做了保护处理。有没有可能某个特征的方差确实是0如果原始数据里某一列全是同一个常数就会有这种情况。我的训练数据里出现过类似问题某个监测指标在多个小时里没有任何变化如果不保护归一化后的结果就是nan后续所有梯度都会变成nan排查起来很痛苦。4.3 梯度下降训练主代码接下来是我自己写的梯度下降实现不调用sklearn。def train_linear_regression(X_train, y_train, learning_rate0.1, epochs2000): num_samples, num_features X_train.shape w np.zeros(num_features) # 初始化权重为0 loss_history [] for epoch in range(epochs): y_pred X_train.dot(w) error y_pred - y_train # 向量长度为num_samples grad (X_train.T.dot(error)) / num_samples w w - learning_rate * grad loss np.mean(error ** 2) / 2 loss_history.append(loss) if epoch % 200 0: print(fepoch {epoch}, loss {loss:.4f}) return w, loss_history这段代码很短但包含了几个关键点权重千万别随机初始化直接全部置零就好。线性回归的损失函数是凸函数从零开始不会陷入局部最优。梯度方向是X^T * error不是error * X^T矩阵维度一定要对。loss除以2这个操作会让打印出来的数字和教科书上的公式一致但梯度里并没有多一个1/2因为误差平方对权重求导之后2会消掉。4.4 在验证集上算RMSE并保存结果训练完成后验证模型效果用的是均方根误差RMSE也就是sqrt(mean((y_true - y_pred)^2))。RMSE的单位和PM2.5相同读起来会很直观。def evaluate_model(w, X_test, y_test): y_pred X_test.dot(w) mse np.mean((y_test - y_pred) ** 2) rmse np.sqrt(mse) # 还可以算平均绝对误差 mae np.mean(np.abs(y_test - y_pred)) return rmse, mae w, loss_history train_linear_regression(X_train, y_train, learning_rate0.1, epochs2000) rmse, mae evaluate_model(w, X_test, y_test) print(f验证集 RMSE {rmse:.4f}) print(f验证集 MAE {mae:.4f})我实际跑出来的验证集RMSE在8到12之间。不同版本的数据集、不同的特征选择结果会有差别。如果你看到RMSE是20多甚至30多先不要怀疑模型公式先检查数据是否错位。我之前有一次RMSE特别高最后发现是滑窗时标签取错了行整个预测序列整体往后平移了一小时导致模型始终在预测“上一小时”的PM2.5误差自然很大。这一步做对以后你可以把权重w保存下来然后对测试集的X_test做同样的归一化再算预测值np.save(linear_regression_weight.npy, w)输出文件可以是CSV里面存储样本编号和预测的PM2.5值。注意输出时不要做反归一化因为我们预测的是PM2.5本身这个值不需要标准化。5. 实测中踩过的坑和几组关键对照实验5.1 训练Loss不降反升问题出在这几个地方我第一版跑出来的loss曲线在几百轮之后突然变成nan当时第一反应是学习率设置过大。降低学习率后变成了“不下降”于是又怀疑是梯度公式写错。排查来排查去最后发现是std[std 0]这一步保护没做好某列方差为0导致归一化后出现无限大值。这件事给我的经验是训练出问题时先检查数据是否包含无穷值或缺失值再检查梯度维度最后才是调学习率。很多新手的顺序是反过来的一看到loss爆炸就忙着把学习率从0.1改到0.00001结果真正的bug根本没解决。我把排查链路整理成下面这个顺序检查X里有没有nan如果有说明归一化或者读取阶段出问题了。检查标签y里有没有特别离谱的负数。PM2.5浓度理论上不会是负值如果出现负值多半是数据列选错了。检查X_train.shape和y_train.shape第一个维度是否一致。不一致的话矩阵乘法的broadcast会悄悄帮你把错误掩盖掉直到算loss时才爆发。打印前几个样本的梯度范数如果数值大于1e6基本可以断定数据没归一化好。5.2 只预测平均值为什么不是放弃思考的理由有一个简单得离谱的baseline不管输入是什么永远预测训练集PM2.5的平均值。这个模型没有使用任何输入特征但它给出的验证集RMSE可以作为下限参考。我计算了一下如果训练集PM2.5均值大约是30左右直接预测均值RMSE可能接近15甚至更高。有了这个baseline我们就能判断线性回归模型是不是真的学到了东西。如果线性回归的RMSE只比“预测均值”好一点点说明模型没有利用好特征。这时候应该去看每个特征的权重值。权重接近0的说明该特征贡献很小权重符号为负的说明该特征与PM2.5负相关。比如风速权重通常是负的风速越大PM2.5越低这是符合直觉的。5.3 全特征和单特征的对照结果我在同一份数据上分别做了两组实验模型特征维度验证集RMSE训练时间baseline预测均值0接近150秒线性回归只用PM2.5九小时910.55秒线性回归全部18特征九小时1628.920秒这个结果可以看出增加特征带来了明显的精度提升但提升幅度不是线性的。说明有些特征之间存在冗余PM10和PM2.5本身高度相关同时放进去会带来共线性问题但总体收益大于共线性的代价。如果你有精力还可以尝试对特征做相关性分析删掉高度相关的指标后RMSE可能不会下降训练更稳定这是后续优化的方向。5.4 关于验证集划分别让邻居来“作弊”空气质量数据有一个很强的特性时间邻近的样本高度相似。如果只用随机划分来评估模型验证集里会混入大量和训练集相隔一小时的样本相当于模型见过几乎一模一样的输入。这样得到的RMSE会偏低但真的在野外遇到新的、未来一小时的数据时表现往往没有验证集上那么好看。我的做法是把最后20%的时间段单独拿出来当验证集模拟模型上线后遇到的未来数据。这样做的风险是如果某个月份空气污染整体偏高验证集和训练集的分布会有差异。但作为作业验证这种“按时间切分”的方式更贴近真实场景也更经得起推敲。6. 做完这道题后的几点延伸思考6.1 线性回归的下一个升级方向是什么线性回归假设输入和输出之间是线性关系但PM2.5的生成、积累、扩散过程本质上是个复杂的非线性过程。做完作业之后很自然的想法是升级到广义线性模型或者给特征加上多项式交叉项。我在实验里试过把PM2.5的历史数据取对数后再做回归发现对高浓度极端值的预测略有改善因为对数压缩了数值范围。如果你有时间可以做这样一个扩展构造“过去9小时PM2.5变化趋势”特征比如当前值减9小时前值、最大值、均值等。这类手工特征对线性模型效果不错因为它把时间变化信息显式提供给模型弥补线性模型不能自动学习时间模式的缺陷。6.2 数据质量的重要性被大多数人低估回过头看这个作业里最花时间的部分不是训练而是处理NR、处理转置后的表头混乱、检查滑窗边界。很多人在分享作业代码时都默认数据是“干净的”但真实世界里的空气质量数据远没有那么规整。传感器可能临时故障、数据传输可能丢失、异常高值可能是附近建筑施工造成的。我个人的体会是如果在数据预处理阶段多花一小时检查训练模型阶段就能少花五小时排错。线性回归模型本身简单可靠不太容易出错出错往往就出在输入数据的形状和语义没有对齐。6.3 最后给你一条实操建议如果你正在做这道作业而且先前完全用调库的方式跑通了我建议你花一个晚上把模型实现换成纯NumPy版本从头写一遍。不为了交作业而是为了亲眼看到权重向量在梯度下降中逐步变化的那个过程。你会发现所谓“机器学习”里那个“学习”的部分不过是一行简单的更新公式在反复执行但就是这行公式能从30万条真实污染记录中总结出一条可复用的预测规律。这种体验是直接调用sklearn拿不到的东西。