李宏毅机器学习作业1:预测PM2.5（含训练数据）

1.要求

给定训练集train.csv，要求根据前9个小时的空气监测情况预测第10个小时的PM2.5含量。

训练集介绍：

CSV文件，包含台湾丰原地区240天的气象观测资料(取每个月前20天的数据做训练集，12月X20天=240天，每月后10天数据用于测试）
每天的监测时间点为0时，1时…到23时，共24个时间节点;
每天的检测指标包括CO、NO、PM2.5、PM10等气体浓度，是否降雨、刮风等气象信息，共计18项；
数据集地址：https://pan.baidu.com/s/1o2Yx42dZBJZFZqCa5y3WzQ，提取码：qgtm。

2. 思路分析

前注：下文中提到的“数据帧”并非指pandas库中的数据结构DataFrame,而是指一个二维的数据包。

2.1 数据预处理

训练集中数据排列形式符合人类观察数据的习惯，但并不能直接拿来喂给模型进行训练，因此需要对数据进行预处理。

浏览数据可知，数据中存在一定量的空数据NR，且多存在于RAINFALL一项。对于空数据，常规的处理方法无非就是删除法和补全法两种。查阅资料后发现，RAINFALL表示当天对应时间点是否降雨，有降雨值为1，无降雨值为NR，类似于布尔变量。因此可以采用补全法处理空数据：将空数据NR全部补为0即可。

根据作业要求可知，需要用到连续9个时间点的气象观测数据，来预测第10个时间点的PM2.5含量。针对每一天来说，其包含的信息维度为(18,24)(18项指标，24个时间节点)。可以将0到8时的数据截取出来，形成一个维度为(18,9)的数据帧，作为训练数据，将9时的PM2.5含量取出来，作为该训练数据对应的label；同理可取1到9时的数据作为训练用的数据帧，10时的PM2.5含量作为label…以此分割，可将每天的信息分割为15个shape为(18,9)的数据帧和与之对应的15个label。

训练集中共包含240天的数据，因此共可获得240X15=3600个数据帧和与之对应的3600个label

# 数据预处理
def dataProcess(df):x_list, y_list = [], []# df替换指定元素，将空数据填充为0df = df.replace(['NR'], [0.0])array = np.array(df).astype(float)#astype就是转换numpy数组的数据类型为float型# 将数据集拆分为多个数据帧for i in range(0, 4320, 18):for j in range(24 - 9):mat = array[i:i + 18, j:j + 9]label = array[i + 9, j + 9]  # 第10行是PM2.5x_list.append(mat)y_list.append(label)x = np.array(x_list)y = np.array(y_list)return x, y, array

x.shape:(3600,18,9)
y.shape:(3600,)

2.2 模型建立

2.2.1 回归模型

采用最普通的线性回归模型，并没有用上训练集中所有的数据，只用到了每个数据帧样本中的9个PM2.5含量值：
y=∑i=08wixi+by=\sum_{i=0}^{8} w_{i} x_{i}+b y=i=0∑8wixi+b
xix_{i}xi为对应数据帧中第i个PM2.5含量，wiw_{i}wi为其对应的权重值，bbb为偏置项，yyy为该数据帧样本的预测结果。

2.2.2 损失函数

用预测值与label之间的平均欧式距离来衡量预测的准确程度，并充当损失函数（这里的损失指的是平均损失；乘1/2是为了在后续求梯度过程中保证梯度项系数为1，方便计算）：
Loss⁡hocl =12num ∑n=0num−1(y^n−yn)2\operatorname{Loss}_{\text {hocl }}=\frac{1}{2 \text { num }} \sum_{n=0}^{ { num }-1}\left(\hat{y}^{n}-y^{n}\right)^{2} Losshocl =2 num 1n=0∑num−1(y^n−yn)2
y^n\hat{y}^{n}y^n为第n个label，yny^{n}yn为第n个数据帧的预测结果，numnumnum为参加训练的数据帧样本个数。

为了防止过拟合，加入正则项：
Loss⁡regularization=12∑i=08wi2\operatorname{Loss}_{\text {regularization}}=\frac{1}{2} \sum_{i=0}^{8} w_{i}^{2} Lossregularization=21i=0∑8wi2

Loss =Losslabel +β⋅Lossregularization =12[1mum⁡∑n=0num−1(y^n−yn)2+β∑i=08wi2]\text {Loss }=L o s s_{\text {label }}+\beta \cdot L o s s_{\text {regularization }}=\frac{1}{2}\left[\frac{1}{\operatorname{mum}} \sum_{n=0}^{n u m-1}\left(\hat{y}^{n}-y^{n}\right)^{2}+\beta \sum_{i=0}^{8} w_{i}^{2}\right] Loss =Losslabel +β⋅Lossregularization =21[mum1n=0∑num−1(y^n−yn)2+βi=0∑8wi2]

β\betaβ为正则项系数。

2.2.3 梯度更新

梯度计算：需明确此时的目标是使Loss最小，而可优化的参数为权重w和偏置值b，因此需要求Loss在w上的偏微分和Loss在b上的偏微分。
∂Loss∂wi=∂Losslabel ∂y∂y∂wi+∂Lossregularization∂wi=1mum∑n=0num−1(y^n−∑i=08wixi−b)⋅(−xi)+β⋅∑i=08wi∂Loss∂b=∂Losslabel ∂y∂y∂b+∂Lossregularization∂b=1mum∑n=0num−1(y^n−∑i=08wixi−b)⋅(−1)\begin{array}{l}{\frac{\partial L o s s}{\partial w_{i}}=\frac{\partial L o s s_{\text {label }}}{\partial y} \frac{\partial y}{\partial w_{i}}+\frac{\partial L o ss_{\text {regularization}}}{\partial w_{i}}=\frac{1}{{mum}} \sum_{n=0}^{num-1}\left(\hat{y}^{n}-\sum_{i=0}^{8} w_{i} x_{i}-b\right) \cdot\left(-x_{i}\right)+\beta \cdot \sum_{i=0}^{8} w_{i}} \\ {\frac{\partial L o s s}{\partial b}=\frac{\partial L o s s_{\text {label }}}{\partial y} \frac{\partial y}{\partial b}+\frac{\partial L o ss_{\text {regularization}}}{\partial b}=\frac{1}{m u m} \sum_{n=0}^{num-1}\left(\hat{y}^{n}-\sum_{i=0}^{8} w_{i} x_{i}-b\right) \cdot(-1)}\end{array} ∂wi∂Loss=∂y∂Losslabel ∂wi∂y+∂wi∂Lossregularization=mum1∑n=0num−1(y^n−∑i=08wixi−b)⋅(−xi)+β⋅∑i=08wi∂b∂Loss=∂y∂Losslabel ∂b∂y+∂b∂Lossregularization=mum1∑n=0num−1(y^n−∑i=08wixi−b)⋅(−1)
计算出梯度后，通过梯度下降法实现参数更新。
wnewi=wi−ηw∂Loss∂wi,bnew=b−ηb∂Loss∂bw_{n e w i}=w_{i}-\eta_{w} \frac{\partial L o ss}{\partial w_{i}}, b_{n e w}=b-\eta_{b} \frac{\partial L o s s}{\partial b} wnewi=wi−ηw∂wi∂Loss,bnew=b−ηb∂b∂Loss
ηw\eta_{w}ηw为权重w更新时的学习率，ηb\eta_{b}ηb为偏置b更新时的学习率。

2.2.3 学习率更新

为了在不影响模型效果的前提下提高学习速度，可以对学习率进行实时更新：即让学习率的值在学习初期较大，之后逐渐减小。这里采用比较经典的adagrad算法来更新学习率。
ηn=ηn−1∑i=1n−1gradi2\eta_{n}=\frac{\eta_{n-1}}{\sqrt{\sum_{i=1}^{n-1} g r a d_{i}^{2}}} ηn=∑i=1n−1gradi2ηn−1
ηn\eta_{n}ηn为更新后的学习率，ηn−1\eta_{n-1}ηn−1为更新前的学习率。∑i=1n−1gradi2\sqrt{\sum_{i=1}^{n-1} g r a d_{i}^{2}}∑i=1n−1gradi2为在此之前所有梯度平方和的二次根。

# 更新参数，训练模型
def train(x_train, y_train, epoch):bias = 0  # 偏置值初始化weights = np.ones(9)  # 权重初始化learning_rate = 1  # 初始学习率reg_rate = 0.001  # 正则项系数bg2_sum = 0  # 用于存放偏置值的梯度平方和wg2_sum = np.zeros(9)  # 用于存放权重的梯度平方和for i in range(epoch):b_g = 0w_g = np.zeros(9)# 在所有数据上计算Loss_label的梯度#Loss在w上的偏微分和Loss在b上的偏微分#参加训练的数据帧样本个数num采用了3200，0~3200是训练集，3200~3600是验证集for j in range(3200):b_g += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-1)for k in range(9):w_g[k] += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-x_train[j, 9, k])# 求平均，除以numb_g /= 3200w_g /= 3200#  加上Loss_regularization在w上的梯度for m in range(9):w_g[m] += reg_rate * weights[m]# adagrad,相当于grad**2的加权和，adagrad算法来更新学习率bg2_sum += b_g ** 2wg2_sum += w_g ** 2# 更新权重和偏置bias -= learning_rate / bg2_sum ** 0.5 * b_gweights -= learning_rate / wg2_sum ** 0.5 * w_g# 每训练200轮，输出一次在训练集上的损失if i % 200 == 0:loss = 0for j in range(3200):loss += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) ** 2print('after {} epochs, the loss on train data is:'.format(i), loss / 3200)return weights, bias

3.完整代码与结果分析

3.1 testPM2.5.py

import pandas as pd
import numpy as np# 数据预处理
def dataProcess(df):x_list, y_list = [], []# df替换指定元素，将空数据填充为0df = df.replace(['NR'], [0.0])array = np.array(df).astype(float)#astype就是转换numpy数组的数据类型为float型# 将数据集拆分为多个数据帧for i in range(0, 4320, 18):for j in range(24 - 9):mat = array[i:i + 18, j:j + 9]label = array[i + 9, j + 9]  # 第10行是PM2.5x_list.append(mat)y_list.append(label)x = np.array(x_list)y = np.array(y_list)return x, y, array# 更新参数，训练模型
def train(x_train, y_train, epoch):bias = 0  # 偏置值初始化weights = np.ones(9)  # 权重初始化learning_rate = 1  # 初始学习率reg_rate = 0.001  # 正则项系数bg2_sum = 0  # 用于存放偏置值的梯度平方和wg2_sum = np.zeros(9)  # 用于存放权重的梯度平方和for i in range(epoch):b_g = 0w_g = np.zeros(9)# 在所有数据上计算Loss_label的梯度#Loss在w上的偏微分和Loss在b上的偏微分#参加训练的数据帧样本个数num采用了3200，0~3200是训练集，3200~3600是验证集for j in range(3200):b_g += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-1)for k in range(9):w_g[k] += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-x_train[j, 9, k])# 求平均，除以numb_g /= 3200w_g /= 3200#  加上Loss_regularization在w上的梯度for m in range(9):w_g[m] += reg_rate * weights[m]# adagrad,相当于grad**2的加权和，adagrad算法来更新学习率bg2_sum += b_g ** 2wg2_sum += w_g ** 2# 更新权重和偏置bias -= learning_rate / bg2_sum ** 0.5 * b_gweights -= learning_rate / wg2_sum ** 0.5 * w_g# 每训练200轮，输出一次在训练集上的损失if i % 200 == 0:loss = 0for j in range(3200):loss += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) ** 2print('after {} epochs, the loss on train data is:'.format(i), loss / 3200)return weights, bias# 验证模型效果
#验证集3400~3600
def validate(x_val, y_val, weights, bias):loss = 0for i in range(400):loss += (y_val[i] - weights.dot(x_val[i, 9, :]) - bias) ** 2return loss / 400def main():# 从csv中读取有用的信息# 由于大家获取数据集的渠道不同，所以数可据集的编码格式能不同# 若读取失败，可在参数栏中加入encoding = 'gb18030'df = pd.read_csv('train.csv', usecols=range(3, 27))x, y, _ = dataProcess(df)# 划分训练集与验证集x_train, y_train = x[0:3200], y[0:3200]x_val, y_val = x[3200:3600], y[3200:3600]epoch = 2000  # 训练轮数# 开始训练w, b = train(x_train, y_train, epoch)# 在验证集上看效果loss = validate(x_val, y_val, w, b)print('The loss on val data is:', loss)if __name__ == '__main__':main()

3.2 结果显示

after 0 epochs, the loss on train data is: 955.3009375
after 200 epochs, the loss on train data is: 49.86823677027294
after 400 epochs, the loss on train data is: 46.20101423801224
after 600 epochs, the loss on train data is: 44.88913061600439
after 800 epochs, the loss on train data is: 44.26903588227097
after 1000 epochs, the loss on train data is: 43.950109190566856
after 1200 epochs, the loss on train data is: 43.78092633274224
after 1400 epochs, the loss on train data is: 43.68982565130423
after 1600 epochs, the loss on train data is: 43.640314303297686
after 1800 epochs, the loss on train data is: 43.61322589236443
The loss on val data is: 40.35422383809947

可以看出，模型在验证集上的损失为40左右，即预测值与label之间的平均差异在6到7之间，由此可见，模型的整体效果还是比较差的。