简介：面向机器学习入门者与数据分析学习者，基于线性回归的PM2.5预测系统完整源码包，围绕空气污染物浓度预测任务，提供从数据读取、预处理到模型训练与结果导出的全流程Python实现。压缩包共19个文件，以12个csv数据文件为主，涵盖训练集、测试集及预测输出；另含5个png结果可视化图、1个npy模型参数文件及1个主程序py脚本，整体大小2.14MB，结构清晰便于直接运行与二次开发。资源中代码对繁体中文big5编码的csv进行了正确处理，并演示了矩阵运算、最小二乘法参数求解、特征拼接等关键步骤，适合课程设计、毕业设计或竞赛练手。目前已有593人浏览学习，可作为理解线性回归实际应用的直观参考，帮助读者快速复现实验并对比预测效果。
1. 用线性回归预测 PM2.5，先分清"能预测"和"预测得准"
PM2.5 浓度预测在环境监测和健康预警场景里一直有真实需求，而它恰恰是线性回归算法教学案例里最"贴地"的一个。PM2.5 与气象变量（温度、气压、湿度、风速）之间存在可解释的线性相关趋势，用线性回归可以在几分钟内得到一条可用的预测基线。这不是玄学，而是能把特征选择、数据切分、残差诊断、正则化这些环节全部落在真实数据上验证。这里要展开的是一套 Python 源码级方案：以公开的北京 PM2.5 数据集为输入，用 scikit-learn 的 LinearRegression 建模，把评估、调优、模型持久化串成可离线运行的预测系统。它适合正在做课程设计、刚入坑机器学习回归任务，或者想把"预测系统"从概念落成可运行代码的工程师。
2. 数据先行：PM2.5 预测的数据集选择与预处理
2.1 公开数据集怎么找：UCI 北京 PM2.5 数据集
做 PM2.5 预测，第一步不是写模型，而是拿到一份字段完整的观测数据。常见做法是使用 UCI Machine Learning Repository 上的 Beijing PM2.5 Data Set，它记录了 2010 年 1 月 1 日到 2014 年 12 月 31 日北京某监测站点逐小时的气象与污染物浓度。数据总量约 4.3 万条，覆盖四年完整季节周期，对线性回归建模来说，这个数据量既不会让训练发散，也不会小到过拟合不可控。
| 字段 | 含义 | 类型 | 在模型里的角色 |
|---|---|---|---|
| year / month / day / hour | 时间戳四要素 | int | 构造时间特征的基础 |
| pm2.5 | PM2.5 浓度（μg/m³） | float | 目标变量 y |
| DEWP | 露点温度（℃） | float | 数值特征 |
| TEMP | 温度（℃） | float | 数值特征 |
| PRES | 气压（hPa） | float | 数值特征 |
| cbwd | 组合风向 | str | 类别特征，需要编码 |
| Iws | 累计风速（m/s） | float | 数值特征 |
| Is / Ir | 积雪 / 降雨时长（小时） | int | 数值特征，稀疏 |
这个数据集的价值在于它同时包含连续特征、类别特征和时间序列特性。其中 cbwd 是四类风向（NE、NW、SE 和代表静风的 cv），cv 类观测往往对应高浓度时段，在建模时要保留这个原始分类，不要合并。需要提醒的是，UCI 版本早期记录里有 pm2.5 的缺失值，这是真实环境监测数据的特点，直接喂给 sklearn 会报错，必须在预处理里消化掉。
2.2 缺失值与时间戳处理：别让 NaN 吃掉模型的 R²
拿到 CSV 之后，我一般先做三件事：把年月日时拼成时间索引、统计缺失分布、查看风向类别的样本量。代码如下：
import pandas as pd import numpy as np df = pd.read_csv('PRSA_data_2010.1.1-2014.12.31.csv') # 拼接标准时间索引，便于排序和构造滞后特征 df['datetime'] = pd.to_datetime(df[['year', 'month', 'day', 'hour']]) df = df.sort_values('datetime').reset_index(drop=True) print(df.isnull().sum()) # 查看各列缺失数量 print(df['cbwd'].value_counts()) # 查看风向类别分布pd.to_datetime 接收的数据帧必须包含 year、month、day、hour 四列，输出一个可直接排序的 datetime 列；sort_values 按时间重排，让后续 shift 操作严格沿时间轴进行。isnull().sum() 能快速暴露缺失集中区，如果 pm2.5 缺失行出现在数据开头，直接用 dropna 删掉即可；如果缺失出现在中间时段，用前向填充 ffill 更符合物理逻辑——PM2.5 浓度属于缓变过程，前一小时的浓度比全天均值更接近当前真实值。
提示：不要用 SimpleImputer 做全局均值填充处理时间序列。均值会把浓度曲线人为拉平，制造出训练集里不存在的时间跳变，最终削弱的恰恰是线性回归对污染过程的还原能力。
对 cbwd 这类类别特征，最简单的方案是 pandas 的 get_dummies 做独热编码。这里要注意 drop_first 参数：四类风向独热后会生成 4 列，但类别间存在共线性，我习惯设置 drop_first=True，只保留 3 个哑变量列，避免特征矩阵里出现完全相关的列导致正规方程不可逆。
2.3 特征工程：滞后特征与周期变量的组织方式
线性回归对特征的组织方式非常敏感。PM2.5 预测里最常用的先验是"过去的浓度影响现在"，因此要提前把浓度序列的滞后项做进特征矩阵：
# 按小时粒度构造多个滞后项，捕捉短时惯性与日周期 for lag in [1, 3, 6, 12, 24]: df[f'pm2.5_lag{lag}'] = df['pm2.5'].shift(lag) df = df.dropna().reset_index(drop=True)shift(lag) 把 pm2.5 整列下移 lag 行，使第 t 行的滞后特征恰好等于 t-lag 时刻的真实浓度。选择 1、3、6、12、24 这几个步长，对应的分别是小时级惯性、半日趋势和 24 小时日循环；滞后 1 小时的系数通常会显著大于其他滞后项，这说明 PM2.5 浓度具有很强的短时记忆。dropna 必须放在最后，否则前 24 行因 shift 产生的 NaN 会进入训练集，连带污染损失函数的计算。
时间字段 hour 不能直接作为数值特征喂给模型。第 23 时与第 0 时数值上相差 23，物理上只相隔 1 小时，线性回归会把这种虚假的数值距离当成真实规律。改进做法是周期编码：
# 周期编码，把小时映射到单位圆上的两点，保留周期性 df['hour_sin'] = np.sin(2 * np.pi * df['hour'] / 24) df['hour_cos'] = np.cos(2 * np.pi * df['hour'] / 24) df['weekend'] = (df['datetime'].dt.weekday >= 5).astype(int)sin/cos 编码把 hour 转换成单位圆上的两个坐标，模型可以学到"午夜到清晨浓度爬升、午后随扩散下降"这类以 24 小时为周期的模式。hour_sin 和 hour_cos 需要同时保留，只用其中一个会把早晚高峰的信息折叠到一起。weekend 是简单的星期特征，因为周末的交通排放与工作日差异明显，这个特征在包含移动源污染的模型里通常有正向贡献。
到这里，特征矩阵已经包含数值特征、周期特征、哑变量和滞后特征，可以作为第 3 章建模的输入。
3. 用 sklearn 跑通线性回归：最小代码与必调参数
3.1 训练集 / 测试集划分：时序数据不能随机打乱
PM2.5 数据是典型的时间序列。如果调用 train_test_split 默认的随机切分，训练集和测试集会混入同一个时间段样本，模型相当于提前看到了未来，测试集上的 R² 虚高但没有部署参考价值。正确做法是按时间顺序切分：
feature_cols = ['DEWP', 'TEMP', 'PRES', 'Iws', 'Is', 'Ir', 'pm2.5_lag1', 'pm2.5_lag3', 'pm2.5_lag6', 'pm2.5_lag12', 'pm2.5_lag24', 'hour_sin', 'hour_cos', 'weekend'] X = df[feature_cols].values y = df['pm2.5'].values cut = int(len(X) * 0.8) # 前 80% 时间段训练，后 20% 时间段验证 X_train, X_test = X[:cut], X[cut:] y_train, y_test = y[:cut], y[cut:]这里刻意不用 sklearn 的 train_test_split，避免有人顺手把 shuffle=True 也带进来。按 cut 切分后，训练集全部位于测试集之前，这样测试集评估的是"模型对未知未来的预测能力"，和真实部署场景保持一致。切分比例 8:2 在这个数据量下比较合理，训练集约 3.4 万条，足够稳定估计 14 个特征的回归系数。
3.2 LinearRegression 的核心参数与拟合逻辑
scikit-learn 的 LinearRegression 是机器学习回归任务里接口最精简的回归器之一，拟合与评估的最小代码可以压缩到一个脚本里：
from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score # 默认 fit_intercept=True，模型带偏置项 model = LinearRegression(fit_intercept=True, positive=False) model.fit(X_train, y_train) y_pred = model.predict(X_test) print('R2:', r2_score(y_test, y_pred)) print('RMSE:', np.sqrt(mean_squared_error(y_test, y_pred))) print('MAE:', mean_absolute_error(y_test, y_pred))fit_intercept 控制是否拟合偏置项 b，在 PM2.5 场景里必须保留 True，因为浓度存在一个非零的背景水平，不拟合截距会把这部分压力强行压到某个特征系数上，扭曲系数含义。positive=False 表示允许系数符号为负，这个默认值不需要改——风速对浓度的作用在静稳天气和强风天气下表现不同，强制非负反而违背真实关系。
训练过程实际求解的是线性回归的正规方程闭式解，数据量约 3.4 万行、特征维度 14，闭式解的矩阵求逆耗时在毫秒级，完全不需要像神经网络那样迭代调学习率。这也是线性回归作为 PM2.5 基线模型的一大优势：可以在几秒内完成全流程验证，把精力留给特征工程和残差诊断。
3.3 评估指标与特征标准化经验
三分钟看评估结果。R²、RMSE、MAE 三个指标最好一起看，它们的含义和行为差异在 PM2.5 场景里非常明显。
| 指标 | 计算公式 | 对错误的敏感度 | PM2.5 场景解读 |
|---|---|---|---|
| R² | 1 - SS_res/SS_tot | 反映整体解释力 | 0.85+ 说明特征对浓度变化的解释力较好，0.6 说明存在明显的模式没学到 |
| RMSE | sqrt(mean((y-y_pred)²)) | 对离群点敏感，平方放大 | 污染爆表时段预测偏差会被放大，适合做安全预警 |
| MAE | mean( | y-y_pred | ) |
特征标准化这个变量，在 LinearRegression 场景下常被误解。如果模型是"裸的"LinearRegression，标准化完全不影响预测结果，因为回归系数会自动缩放以抵消特征量纲差异，预测值不变。但一旦引入 Ridge 或 Lasso，正则化项对系数整体大小施加惩罚，量纲大的特征（如 PRES 气压值上千）会被无差别压低，这时必须先标准化再训练，才能保证每个特征受到的惩罚强度一致。
from sklearn.preprocessing import StandardScaler from sklearn.pipeline import make_pipeline from sklearn.linear_model import Ridge # 标准化 + 岭回归组合，避免量纲影响正则化 model_ridge = make_pipeline( StandardScaler(), Ridge(alpha=1.0) ) model_ridge.fit(X_train, y_train)make_pipeline 的作用是把标准化和回归器串成一个整体对象，fit 时先做标准化再做岭回归，predict 时会自动重复同样的变换，避免在测试集上单独调用 scaler.transform 时忘记使用训练集得到的均值和方差。alpha=1.0 是最常见的起调值，之后可以按 0.1、1、10 的对数刻度搜索，观察验证集 R² 的峰值落在哪个区间。
4. 线性回归模型诊断与特征优化：让 R² 从 0.6 走到 0.9
4.1 残差分析：线性回归失效的信号
基线模型跑完，第一个动作永远是看残差。残差是 y_test 减去 y_pred，如果残差对预测值呈现出漏斗形或弯曲的分布，说明模型存在系统性的低估或高估，公式层面没问题也掩盖不了。
import matplotlib.pyplot as plt residual = y_test - y_pred # 残差图：横轴预测值，纵轴残差，观察是否随机散落 plt.scatter(y_pred, residual, s=4, alpha=0.4) plt.axhline(y=0, color='red', linestyle='--') plt.xlabel('Predicted PM2.5 (μg/m³)') plt.ylabel('Residual (μg/m³)') plt.savefig('residual.png', dpi=120)PM2.5 数据残差图几乎必然会展示出两种特征：一是预测值超过 300 时残差发散，真实高浓度时段被明显低估；二是低浓度段残差贴近零轴，说明模型对清洁天气拟合得不错。第一个现象的物理背景是：重污染通常由静稳天气加高湿度共同触发，线性回归只能学到这些变量的平均效应，无法单独刻画极端组合条件。这类问题靠换模型不如靠补特征和调滞后窗口。
4.2 用系数可解释性做特征筛选
线性回归的可解释性在预测类任务里是稀缺品。训练完成后，把系数输出成表格，按绝对值排序，能直接看出模型在依赖哪几个特征：
coef_df = pd.DataFrame({ 'feature': feature_cols, 'coef': model.coef_ }).sort_values('coef', key=abs, ascending=False) print(coef_df.head(10))实际训练结果里，pm2.5_lag1 的系数通常远高于其他特征，这符合 PM2.5 的短时记忆效应；风速 Iws 的系数显著为负，也符合"风大浓度低"的物理直觉。如果某个特征系数绝对值非常小，且训练集和验证集上的表现没有明显差异，就可以把它从 feature_cols 中移除。特征数量从 14 降到 10 后，模型复杂度下降，残差分布往往更稳定。
4.3 多项式扩展：线性回归的另一种升级路径
当四个气象变量与浓度之间不只是简单线性关系时，可以引入交互项和平方项，让线性回归获得非线性表达能力。这一步常见做法是直接上 PolynomialFeatures。
from sklearn.preprocessing import PolynomialFeatures # degree=2 生成一次项、平方项与两两交互项 poly = PolynomialFeatures(degree=2, include_bias=False, interaction_only=False) X_poly = poly.fit_transform(X_train) X_test_poly = poly.transform(X_test) model_poly = LinearRegression() model_poly.fit(X_poly, y_train) y_pred_poly = model_poly.predict(X_test_poly)degree=2 时，特征维度会从 14 膨胀到约 100 个，原因是公式 n_features * (n_features + 1) / 2 涵盖了所有平方项和交叉项。维度大幅升高后，LinearRegression 的方差会上升，容易在测试集上出现 R² 不升反降。这时的正确组合是第 3 章提到的 Ridge，L2 正则化能约束交互项系数的范数，让模型保留多项式表达能力的同时不剧烈过拟合。interaction_only=True 则只生成交互项，不给平方项，适合已有明确先验"温度与湿度的乘积影响浓度"的场景。
一个必须警惕的坑：PolynomialFeatures 必须先在训练集上 fit，再对测试集单独 transform，不能对全量数据做扩展后再切分。因为拟合过程会记录每列的均值、方差和特征名称，测试集信息一旦混入 fit 阶段，评估结果就会虚高。这一点和 StandardScaler 的使用规则完全一致。
4.4 用 TimeSeriesSplit 验证跨时段稳定性
单次 8:2 切分只能说明模型在某一特定未来时段的表现。要验证模型在跨季节、跨年份是否稳定，应该用时间序列交叉验证。sklearn 的 TimeSeriesSplit 专门解决这个问题：
from sklearn.model_selection import TimeSeriesSplit, cross_val_score # 5 折滚动验证：每次用前 k-1 折训练，第 k 折验证 tscv = TimeSeriesSplit(n_splits=5) scores = cross_val_score(model, X, y, cv=tscv, scoring='r2') print(scores, scores.mean())TimeSeriesSplit 的折切分永远保持训练集在验证集之前，n_splits=5 会把全量数据切成 5 个连续时段，依次做 5 次训练与验证。输出结果如果显示出前几折 R² 高、后几折 R² 低，说明模型对近期数据的拟合能力在衰退，原因是污染物组成和气象条件在年份间有缓慢漂移。反过来的趋势说明数据尾部存在极端污染事件，模型预测能力被某个时段拉低。
5. 交付离线预测系统：模型持久化与定时预测脚本
5.1 joblib 保存与加载模型的正确姿势
模型调好之后，训练脚本和预测脚本要分离。sklearn 官方推荐用 joblib 持久化模型对象，它比 pickle 对大对象序列化更高效：
import joblib # 保存训练好的模型到本地文件 joblib.dump(model, 'pm25_lr_model.joblib') # 加载并直接用于预测 loaded_model = joblib.load('pm25_lr_model.joblib')joblib.dump 的第一个参数是模型对象，第二个参数是保存路径。加载后的模型不需要重新 fit，predict 方法与原先完全一致。需要提醒的是，保存时顺手把 feature_cols 也序列化成一个 JSON 文件或者 joblib 文件，加载预测时通过字段对比校验输入顺序，能避免因后续新增特征而导致的"预测时列错位"问题。
5.2 一个可直接排期的预测脚本骨架
预测系统的最终形态，通常是一个每小时被 cron 调起的脚本。它读取最近一个小时的观测值，构造特征，输出下一小时浓度预测，并把结果追加写入日志：
import pandas as pd import numpy as np import joblib # 加载模型，全局只加载一次 model = joblib.load('pm25_lr_model.joblib') def build_features(obs): # 输入为一条最新观测记录，输出与训练特征顺序一致 return [obs['DEWP'], obs['TEMP'], obs['PRES'], obs['Iws'], obs['Is'], obs['Ir'], obs['pm2.5_lag1'], obs['pm2.5_lag3'], obs['pm2.5_lag6'], obs['pm2.5_lag12'], obs['pm2.5_lag24'], np.sin(2 * np.pi * obs['hour'] / 24), np.cos(2 * np.pi * obs['hour'] / 24), int(obs['weekday'] >= 5)] latest = pd.read_csv('latest_weather.csv').iloc[-1] features = build_features(latest) pred = model.predict([features])[0] # 追加写入预测结果，供后续自检 with open('predictions.log', 'a') as f: f.write(f"{latest['datetime']},{pred:.1f}\n")build_features 函数里的字段顺序必须与训练时的 feature_cols 完全一致，这是预测脚本最容易出错的地方。latest_weather.csv 由上游的数据采集服务写入，脚本只读最后一行，不关心采集细节，这样的设计把数据流向拆成"采集 → 预测 → 落盘"三段，任何一段都可以独立替换。追加写日志的方式简单但可靠，之后可以用同一组历史数据画时间序列曲线。
注意：cron 每小时调用脚本时，需要保证滞后特征使用的历史浓度是已观测到的真实值，而不是上一轮的预测值。使用预测值做滞后输入，误差会在小时尺度上发生累积，也就是常说的"反馈误差漂移"。
5.3 预测结果自检：用最后 7 天滑动 MAE 判断模型漂移
系统部署后，我保留一个固定位置的最近 7 天窗口，每天把模型输出的预测值与真实监测值做一次对比，算滑动平均绝对误差，超过阈值就告警。原因是环境监测传感器会缓慢漂移，气象输入分布也会随季节改变，模型第一天表现好不代表第二个月仍然稳定。
df_eval = pd.read_csv('predictions.log', names=['datetime', 'pred']) df_real = df_recent[['datetime', 'actual']] merged = pd.merge(df_eval, df_real, on='datetime', how='inner') merged['err'] = (merged['pred'] - merged['actual']).abs() # 输出最近 7 天（168 小时）的平均绝对误差 print(merged.tail(7 * 24)['err'].mean())当结果高于 35 μg/m³，先把最近几天大风、静稳等天气事件标出来核对；如果误差集中出现在固定时段，优先怀疑数据采集端的时间戳对齐问题，而不是模型本身。用这样一个轻量自检脚本，离线预测系统就具备了基本的可观测性，后续再往里面加新的污染源特征，也能用同样的滑动窗口做 A/B 对比，评估新特征的真实增益。
本文还有配套的精品资源，点击获取
