这个题目我盯了好一会儿原因很简单这是我当年认认真真从零手写、一行库都没调过的机器学习作业。很多人觉得课程作业嘛调个包跑通就行但那次我硬是拿numpy把线性回归、梯度下降、特征工程全部手撸了一遍。结果就是后面再遇到各种看似高深的模型我一点都不慌因为最底层的计算逻辑已经被我亲手拆开揉碎过一遍了。这篇就完整复盘一下这个“预测PM2.5”的作业题怎么拆、数据怎么处理、模型怎么选、代码怎么写、坑怎么踩最后再把一些可以迁移到其他项目的经验一并交代清楚。无论你是刚开始学机器学习的学生还是想补基础底子的转行者这篇应该都能帮你省下不少摸索的时间。1. 题目拆解与方案选型为什么是线性回归为什么非要手写1.1 题目的真实意图练的不是预测是原理先说作业本身。数据集是某高校机器学习课程提供的内容是某地区一整年的气象观测数据每行记录了一个观测站某个小时的18项指标包括PM2.5浓度、PM10浓度、NO2、CO、O3、风速、雨量、日照等。任务是给定前9个小时的观测数据预测第10个小时的PM2.5值。看起来是一个时间序列预测问题但它其实是课程用来练线性回归的。为什么选线性回归因为它是最基础、最直观、也最能暴露你“到底懂不懂机器学习”的模型。逻辑回归适合分类树模型适合非线性但线性回归把“特征加权求和加偏置再最小化误差”这个过程展现得干干净净。梯度下降、损失函数、过拟合、正则化这些核心概念都能在线性回归上找到最直观的落脚点。我当时还纠结过要不要直接上LSTM或者XGBoost后来想想不对作业要的是让你理解模型参数是怎么被训练出来的不是让你炫技。用10行代码调一个xgboost跑出好看的分数对你理解机器学习帮助为零。手写线性回归虽然笨但它逼着你去面对每个细节梯度怎么算、学习率该设多少、特征要不要归一化、数据怎么切才不会偷看未来信息。1.2 手写代码和调库的本质差别很多同学觉得作业用sklearn的LinearRegression跑一下再report一个RMSE就完事了。不是说不行而是这样你永远只能在“调包侠”的水平线上打转。手写代码强迫你回答几个关键问题梯度下降的更新公式为什么是w - learning_rate * gradient而不是w learning_rate * gradient损失函数对w求导时那个转置到底怎么来的为什么要转置数据归一化后预测结果怎么还原成原始量纲正则化项加在损失函数里梯度的表达式会发生什么变化这些问题你在sklearn里永远碰不到因为fit方法封装得干干净净。但真实工作中你迟早会遇到“现成库跑不了”的场景那时候你能依赖的只有原理本身。我自己的体会是手写一遍等于把“黑盒模型”变成“灰盒模型”后面学神经网络、学深层模型脑子里始终有那张矩阵乘法的图。1.3 方案选型的思考链路我最终确定的方案是纯numpy手写线性回归用批量梯度下降Batch Gradient Descent训练搭配L2正则化防止过拟合特征用StandardScaler做标准化评价指标用RMSE和MAE。选这个组合的理由相当朴素数据量不大不到4000条样本18个特征批量梯度下降完全跑得动没必要上SGD或Mini-batch增加复杂度。加上L2正则化是因为特征经过多项式扩展后维度会膨胀不加正则很容易在训练集上自嗨。不用正规方程Normal Equation因为梯度下降才是课程的核心考点而且正规方程在特征维度高的时候计算逆矩阵代价很大不利于后续扩展。这个选择在当时看起来偏“笨”但后来我意识到批量梯度下降其实帮你自动完成了对学习率的敏感性测试这对理解优化算法非常有帮助。2. 数据预处理与特征工程把原始表格变成能喂给模型的模样2.1 原始数据的真实长相坦白讲第一次拿到数据的时候我是有点懵的。CSV文件里根本不是标准的一行一个样本而是按“日期小时”的行组织方式。每一天有24行每行有18个特征列然后后面还有第19列是当小时的PM2.5测量值。这意味着如果整表直接喂给线性回归模型会把“同一日期的第14行”和“第15行”当成完全独立的样本丢失了时间顺序信息。要预测第10小时的PM2.5需要把前9个小时的数据拼接成一个长的特征向量。也就是说样本的构造方式要从“逐行预测”改成“按滑动窗口拼接”。我当时设计了一个窗口函数简单说就是把连续9个小时的18维特征拼接成162维向量然后对应第10个小时的PM2.5作为标签。这样构造出来的数据集每一行才是真正意义上的独立样本。这个思路和现在所有时间序列预测模型的数据预处理方向一致只是当时做起来很原始。2.2 缺失值处理的三个备选方案数据里有不少RAINFALL字段是空值NR表示没有降雨。直接填0其实问题不大因为雨量是连续值没下雨用0代替是合理的。但更规范的方案是看这个特征与目标值的相关性如果当前列缺失占比过高干脆丢弃这一列否则就用该列均值填充。我的做法是对数值型特征列统计缺失比例缺失超过50%的直接删掉。缺失比例较低的列用该列均值填充。RAINFALL字段统一填0因为“NR”在这里语义明确就是无雨。这里有个小坑容易被忽略填充均值时必须使用训练集的均值不能把测试集也统计进去。数据和标签分开前你就把整个数据集的均值算好塞到测试集里那就是标签泄露会让分数虚高但实际部署时完全失效。我自己第一次就没注意后来验证集分数和公有分数差距很大排查半天才发现是这个问题。2.3 标准化为什么线性回归一定要做这一步线性回归的特征量纲差异极大。PM2.5浓度可能只有几十而CO浓度可能只有0.5左右但风速可能是几十气压可能上千。如果不做标准化梯度下降会走得非常痛苦某些维度的梯度大得离谱某些维度梯度小得可怜整个训练过程像瘸腿走路不仅收敛慢还容易震荡不收敛。标准化公式很简单x_scaled (x - mean) / std用numpy实现也就是两行代码的事。但要注意三点要保存训练集的mean和std测试集和真正预测时都用同一组参数来变换不能各自算各自的。预测出标准化后的PM2.5值后要反过来用训练集标签的mean和std还原成真实数值这样才能算RMSE。标准化不仅对输入特征有必要对标签做标准化也常常能提升数值稳定性。我当时把标签也标准化了训练时损失函数的初始值会变得很小收敛速度肉眼可见地加快。2.4 特征扩展为什么把9小时拼一起还不够只把前9小时的18维特征拼成162维线性回归能捕捉的信息其实有限。PM2.5的浓度变化受风力、湿度、前一天浓度共同影响这些关系未必是线性的。所以我又做了两件事把第9小时的PM2.5单独拿出来乘以一个较大的权重因为显然最近一小时的浓度对下一小时有极强的影响。对部分特征尤其是PM2.5、PM10、CO计算二阶交叉特征例如PM2.5和风速的乘积。风速大时污染物扩散快这种交互效应用乘积特征才能表达出来。特征扩展后维度从162变成了接近200配合L2正则化完全能hold住。这里我的体会是特征工程不是越复杂越好但要基于对业务的理解做合理扩展。你可以不懂大气科学但你至少要知道“风速越大污染物浓度越低”这种常识它就是一组有意义的交叉项。3. 核心算法原理手写梯度下降时脑子里的那张图3.1 线性回归的目标函数和损失函数线性回归的核心假设是目标值y是特征向量x的线性组合加上一个偏置b也就是y_pred w^T x b为了表达方便通常把b并进w里在x前面加一列全1这样可以直接写成矩阵相乘。损失函数用均方误差MSE公式是L (1/n) * sum((y_pred - y)^2) lambda * sum(w^2)后面那项是L2正则化项。为什么加它因为在样本量相对特征维度不够大的时候模型很容易记住训练集的噪声导致在验证集上表现糟糕。加了正则项之后模型被强迫往“权重尽量小”的方向走相当于给模型套了缰绳。3.2 梯度推导那个转置到底哪来的手写代码真正考验人的就是梯度计算。对损失函数求w的偏导dL/dw (2/n) * X^T (y_pred - y) 2 * lambda * w这里的关键是理解X^T为什么要乘残差向量。可以这样想X的形状是(n, d)残差是(n, 1)X^T就是(d, n)矩阵乘完之后得到(d, 1)正好和w的形状对得上。如果不小心写成X乘以残差得到的是(n, 1)维度完全错乱程序直接报错。我当时是硬推导了两遍才想明白这个转置。实际写代码时更安全的做法是先算出预测值、残差然后再算梯度而不是直接在纸上推完就往代码里抄。推荐新手先在代码里print出梯度矩阵的shape确认是(d, 1)再做减法。3.3 批量梯度下降的完整更新过程训练过程写成伪代码就是初始化 w 为零向量 repeat: y_pred X w error y_pred - y gradient (X.T error) / n lambda * w w w - learning_rate * gradient until 损失变化小于阈值 或 达到最大迭代次数这里有一个细节梯度里正则化项通常不对偏置b做惩罚所以如果w和b是分开存的话要记得给b的梯度去掉正则项。我为了方便直接把b塞进w的第一列这样所有维度一起更新虽然严格来说不够严谨但实际影响非常小尤其是在样本量不小的情况下。学习率的选择可以先用一组小实验来试分别用0.1、0.01、0.001跑100轮看损失曲线是震荡、平稳下降还是几乎不动。我最后选定的学习率是0.001迭代2000轮损失曲线下降得非常平稳接近收敛时也没有明显震荡。4. 实操过程与关键代码实现从零到能提交的完整方案4.1 数据读入与预处理代码整个流程的第一部分是读数据、造特征、切训练验证集。我用的就是pandasnumpy的组合没有调用任何机器学习库。核心代码大致长这样import numpy as np import pandas as pd data pd.read_csv(train.csv, encodingutf-8) # 原始数据中 RAINFALL 的 NR 需要替换为数值 0 data.replace(NR, 0, inplaceTrue) data data.astype(float) # 按日期分组每天24个小时18个特征 months data.groupby([month, day]) def build_sample(group): # group 是某天24行的数据 arr group.iloc[:, 2:].values # 去掉前两列标识列 window [] for hour in range(9): window.append(arr[hour]) feature np.concatenate(window) label arr[9][-1] # 第10个小时的PM2.5 return feature, label features [] labels [] for _, group in months: f, l build_sample(group) features.append(f) labels.append(l) X np.array(features) y np.array(labels)这段代码有个地方要非常小心按日期分组时如果原始数据里某个日期缺了若干个钟头拼接出来的窗口就会错位。在我自己处理的过程中数据整体比较规整没有出现大问题但严谨起见可以在build_sample里面加一个长度校验如果当天不足24行就直接跳过或者报错。最后再把X,y划分成训练集和验证集比例大概是8:2。切分要打乱顺序吗这里要特别注意。因为样本是按日期顺序构造的直接按顺序切会让训练集全是前半年的数据验证集全是后半年的导致季节分布不均衡。我当时的做法是先随机打乱所有样本的索引再切分。这是完全合理的因为我们的目标是预测“任意时刻的下一小时”不是预测“未来一年”所以样本独立性假设成立。4.2 标准化与特征扩展的落地实现特征扩展的逻辑放在标准化之后做还是之前做我建议先标准化再做交叉特征。原因是如果不先标准化降水量几千的数值和PM2.5几十的数值乘起来会得到一个天文数字严重拖慢梯度下降。标准化之后再相乘数值范围就合理多了。mean X_train.mean(axis0) std X_train.std(axis0) std[std 0] 1 # 防止某列方差为0导致除零 X_train_scaled (X_train - mean) / std X_val_scaled (X_val - mean) / std # 添加交互特征 X_train_inter np.hstack([ X_train_scaled, X_train_scaled[:, [8]] * X_train_scaled[:, [6]], # PM2.5与风速 ]) X_val_inter np.hstack([ X_val_scaled, X_val_scaled[:, [8]] * X_val_scaled[:, [6]], ])这里的索引是我临时编的实际要看你数据里第几列是PM2.5、第几列是风速别直接照抄。做特征扩展前建议把列名打印出来一一核对避免乘了一个莫名其妙的特征组合。交叉特征不需要做太多选两三个业务上说得通的就行。加多了不仅训练变慢还容易引入噪声。4.3 手写线性回归模型主体这是整个作业的核心也是我最有感情的一段代码。所有逻辑都在一个类里清晰又方便测试class LinearRegressionManual: def __init__(self, lr0.001, n_iters2000, reg_lambda0.1): self.lr lr self.n_iters n_iters self.reg_lambda reg_lambda self.w None def fit(self, X, y): n_samples, n_features X.shape self.w np.zeros(n_features) X_b np.hstack([np.ones((n_samples, 1)), X]) # 增加一列1 self.w np.zeros(n_features 1) for i in range(self.n_iters): y_pred X_b self.w error y_pred - y gradient (X_b.T error) / n_samples # 不对偏置项加正则化 gradient[1:] self.reg_lambda * self.w[1:] self.w - self.lr * gradient if i % 200 0: loss np.mean(error ** 2) self.reg_lambda * np.sum(self.w[1:] ** 2) print(fiter {i}, loss {loss:.4f}) def predict(self, X): X_b np.hstack([np.ones((X.shape[0], 1)), X]) return X_b self.w这段代码你拿去改改学习率和迭代次数就能用。有几点心得损失打印不需要每轮都print否则控制台刷屏不说还拖慢训练速度。每200轮打印一次刚好。梯度更新的时候正则项只加在gradient[1:]上这和前面说的偏置不惩罚一致。虽然影响很小但这是思路严谨性的体现。如果loss出现NaN优先检查特征里是否有无穷大值或者学习率过大导致数值发散。4.4 验证与评价指标的实现训练完之后要分别计算训练集和验证集的RMSE和MAE。这俩指标的意义不同RMSE对大误差非常敏感一个极端预测能拖垮整体分数MAE更鲁棒反映平均偏差水平。我建议同时看两个如果RMSE明显大于MAE说明模型对异常值敏感可以考虑对标签做截断或者换更稳健的损失函数。y_pred_train model.predict(X_train_inter) y_pred_val model.predict(X_val_inter) # 还原标签标准化 y_pred_train_real y_pred_train * y_std y_mean y_pred_val_real y_pred_val * y_std y_mean rmse_train np.sqrt(np.mean((y_pred_train_real - y_train_real) ** 2)) rmse_val np.sqrt(np.mean((y_pred_val_real - y_val_real) ** 2)) print(fTrain RMSE: {rmse_train:.3f}) print(fVal RMSE: {rmse_val:.3f})如果训练集RMSE很低但验证集RMSE高就是过拟合了。我当时训练集RMSE大概在12左右验证集在16左右差距不算大说明特征扩展配合正则化把握得还可以。如果差距超过30%就该减少特征维度或增大正则化系数。5. 常见问题与排查技巧实录那些年踩过的坑5.1 损失函数爆炸直接变NaN怎么办这个问题几乎每个手写梯度下降的人都遇到过。表现是训练到第几十轮时loss突然变成nan然后永远回不来。原因通常是两种学习率太大导致更新步长越过最优解或者数据里有异常值导致梯度计算数值溢出。我当时排查的顺序是这样一个问题一个问题的验证先打印前10轮的loss看看是初始就nan还是训练到后面才nan。如果是初始就nan检查数据预处理是否有inf或空值没替换干净。如果是后面才nan把学习率除以10重跑一遍如果不再nan说明学习率过高。再把损失函数里的正则项系数从0.1调成0.01排除正则项梯度过大的可能。实际结果我那次是学习率设了0.01对归一化后的数据来说偏高。一路降到0.001才稳定。经验是标准化的数据学习率从0.001起步观察曲线后再考虑增大。不要一上来就设0.1除非你想看烟花。5.2 验证集分数虚高但提交分数崩塌数据泄露这个问题比NaN更阴险因为它表面上损失曲线完美结果也漂亮但一上真正的测试集就露馅。我当时犯的错误就是在构造样本时把第10个小时的数据混进了前9小时的特征里。怎么发生的我在分组处理时为了图方便把当天的18列全部塞进了特征矩阵然后标签用的也是同一行的PM2.5值。这样模型看到的不只是前9小时而是第10小时的PM2.5直接作为特征通向输出。训练时RMSE能低到逼近0但测试集上直接原形毕露。排查方法其实很简单你把训练好的模型权重打出来如果发现某个特征的权重异常大几百上千说明该特征和标签几乎共线八成是数据泄露了。5.3 训练慢得像蜗牛是代码问题还是数据问题还有一次我发现训练速度极慢每轮都要跑很久。后来发现是我的梯度计算用了循环遍历每一条样本而不是矩阵运算。numpy的矩阵乘法底层是经过优化的C代码手写for循环跑几千条样本每轮都要花掉几百毫秒。改成X_b.T error之后训练直接从几分钟降到了十几秒。经验教训是能用矩阵运算就绝不用循环。这也是手写代码过程中最大的性能提升点。你要是把循环版交上去助教大概率会说“功能正确但效率堪忧”。5.4 常见问题速查表现象可能原因排查/解决方法loss初始就为NaN数据有空值/inf未处理检查预处理替换空值和无穷值loss训练中爆炸学习率过大调小学习率或对数据标准化训练慢用了for循环而非矩阵运算检查梯度计算改为矩阵乘法训练集完美但验证集差特征泄露或严重过拟合检查特征是否混入标签增加正则化验证集还行但提交崩数据切分没打乱类别不均衡打乱数据再切分确认同分布权重数值异常大特征未标准化或数据泄露标准化特征检查特征与标签关系预测结果全一样模型退化或偏置太大检查特征方差确认梯度是否更新6. 从作业到实战这几次手写到底给我留下了什么6.1 直接能迁移到工程里的几个习惯手写完成这个作业后我最大的收获不是“我会调线性回归了”而是建立了几个条件反射式的习惯第一个习惯是拿到数据先看形状和缺失值不急着建模。任何一个工业级项目数据的脏乱程度都远超课程作业。先打印info、describe确认列数、行数、缺失值占比再动手预处理能省掉后面90%的排查时间。第二个习惯是每次训练都要划分留出验证集并且从头到尾只用训练集的统计信息来预处理数据。无论是均值、标准差还是缺失值填充测试集和验证集都要复用训练集的那组参数。这一步做对了你的模型泛化能力才谈得上可信。第三个习惯是监控训练指标的曲线而不是只看最终分数。损失函数是否平滑下降、是否震荡、是否突然跳变这些曲线里写满了训练过程的健康状况。很多时间序列预测项目里只盯着一个最终RMSE很容易被数字欺骗曲线才是诊断工具。6.2 如果重新做一次我会改进哪里要是让我现在重新做这个作业我会在特征和模型两个方向上作点升级特征方面加上滞后N小时的PM2.5均值、滚动平均、小时周期编码比如把小时数转成sin和cos分量来捕捉昼夜节律。这些对空气质量预测确实有真实增益而且实现起来并不复杂。模型方面可以试试加一层特征交叉后的岭回归或者用交叉验证来调正则化系数而不是靠手试。另外一个值得尝试的方向是把9小时窗口改成可变长度窗口对比不同窗口长度下模型的表现。窗口太短捕捉不到趋势太长则引入冗余信息。选择一个在验证集上表现最好的窗口长度本身就是一种特征选择。6.3 给后来者的几句实在话如果你也要做类似的课程作业我不建议你直接复制粘贴任何人的代码包括我这篇里的示例。你可以参考思路、看懂结构但一定要亲手敲一遍。因为手写代码的过程中你最缺的不是“会调用某个函数”而是“知道它在干什么”。同样一段归一化代码自己敲一遍才会知道std为0的列会导致除零同样一段矩阵梯度自己推导一遍才清楚转置从哪里来。还有个非常实用的建议每实现一个功能就写一小段测试来验证它。比如标准化之后检查每列均值是否接近0方差是否接近1训练完之后手动计算几条样本的预测值和模型输出对比确保没算错。这些看似繁琐的检查能帮你把错误拦截在早期。最后别害怕调参。学习率、正则化系数、迭代次数这些参数没有绝对的标准答案。用手试用曲线验证你会慢慢培养出一种直觉loss下降太快可能是学习率过大下降太慢要加大学习率震荡则说明步子迈大了。这种直觉不是看书看出来的是亲手调出来的。回头看不就是一次普普通通的课程作业吗但它确确实实帮我打下了机器学习的基础。希望这篇复盘也能帮你在同样的路上走得顺一点。