[机器学习-回归实例] 李宏毅机器学习作业一：PM2.5 预测实例

[机器学习-回归实例] 李宏毅机器学习作业一：PM2.5 预测实例

简介

  这是一个回归预测的实例，数据集来自李宏毅的《机器学习》课程的第一次作业，可以用来训练一个线性回归模型或者其他的回归模型，动手实践学习到的理论知识，在实践中去发现问题。这篇文章主要包括数据集的解释、数据的预处理过程及线性回归模型的代码实现这三个部分。由于这个数据集并不能直接使用，还需进行一些重构，为了对初学者更友好一些，我会提供重构完的可以读取后就直接使用的数据集，我也会提供 $\mathrm{M}\mathrm{a}\mathrm{t}\mathrm{l}\mathrm{a}\mathrm{b}$ 可直接使用的 $.\mathrm{m}\mathrm{a}\mathrm{t}$ 文件供使用 $\mathrm{M}\mathrm{a}\mathrm{t}\mathrm{l}\mathrm{a}\mathrm{b}$ 的同学使用，不过还是希望同学们能够使用原始数据自己动手重构，重构的代码也会附在文中。

数据集解释

  这个数据集是一个城市 $2014$ 年一整年的空气质量监测的实测数据，包含了 $\mathrm{A}\mathrm{M}\mathrm{B}\_\mathrm{T}\mathrm{E}\mathrm{M}\mathrm{P}、\mathrm{C}\mathrm{H}4$ 等十八种关于空气质量的参数，每个小时采集一次。在 $\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}.\mathrm{c}\mathrm{s}\mathrm{v}$ 训练集数据集中包含了每个月前 $20$ 天的每个小时采集的数据，一共有 $12\times 20\times 24$ 个小时的监测数据。 下图截取了两天的前 $10$ 个小时的原始数据的图，可以看到从第 $4$ 列开始每一列是 $1$ 个小时采集的数据，一共有 $24$ 列。从第二行开始每 $18$ 行是一天的 $18$ 种监测数据。

  我们要做的任务回归预测任务就是依据前 $9$ 个小时的 $18$ 种监测的参数预测第 $10$ 个小时的 $\mathrm{P}\mathrm{M}2.5$ 的取值。

数据集重构

  直接用这个数据集显然是不容易进行训练的，因为属性值 ${\mathbf{X}}_{train}$ 和标签 ${\mathbf{y}}_{train}$ 并没有分开，而且在 ${\mathbf{X}}_{train}$ 中应该每一行是一条数据，所以我们数据重构的任务就是从这个原始的数据集 $\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}.\mathrm{c}\mathrm{s}\mathrm{v}$ 来构建我们方便使用的 ${\mathbf{X}}_{train}$ 和 ${\mathbf{y}}_{train}$ 。

  首先来计算一下我们总共可以构建出多少条数据，由于构建数据集的数据在时间上是要连续的，而每个月只有前 $20$ 天的数据，所以我们没法跨月构建数据集，在任意一个月内，我们有 $20(天)\times 24(小时)$ 共 $480$ 小时的数据，所以每 $9$ 小时可以构建 ${\mathbf{X}}_{train}$ 上的一条数据，第 $10$ 个小时的 $\mathrm{P}\mathrm{M}2.5$ 可以作为这一条数据的标签。所以每个月可以有 $471$ 条数据，$12$ 个月就有 $5652$ 条数据。接下来是构建数据集的过程（以下代码是使用 Python 编写实现的，关于环境的配置等如有问题可在评论中联系我，我可以写个配置环境的教程）：

Step1：读取数据集

  调用 $\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{s}$ 包和 $\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{p}\mathrm{y}$ 包，使用 pd.read_csv()函数读取 $\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}.\mathrm{c}\mathrm{s}\mathrm{v}$ 文件。

#读取数据集
import pandas as pd   
import numpy as np

data = pd.read_csv('./train.csv', encoding = 'big5')  #train.csv文件要与python代码文件(.py或.ipynb)在同一个目录下才行

  去除表格中的前三列（第一行自动被去除了，不用主动去除），并将表格中 $\mathrm{R}\mathrm{A}\mathrm{I}\mathrm{N}\mathrm{F}\mathrm{U}\mathrm{L}\mathrm{L}$ 的 $\mathrm{N}\mathrm{R}$ 变成0（因为线性回归不能处理文字数据）。最后使用 to_numpy()函数把数据格式变为 $\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{p}\mathrm{y}$ 的格式。

data = data.iloc[:, 3:]  #去除前三列
data[data == 'NR'] = 0  #将NR全变成0
raw_data = data.to_numpy()  #将数据转换成numpy数据

Step2：按月重构数据

  因为一个月内 $20$ 天的数据在时间上是连续的，前一天的第 $23$ 小时后面的一小时的数据就是后面一天的第 $0$ 小时的数据，而在原来的表格中，每天的数据是在不同列的，所以现在所要做的就是将一整个月的数据连成“一行”（其实是 $18$ 行，因为有 $18$ 种监测的数据），示意图如下所示：

  用一个字典来装所有月的数据，一个列表装具体某一个月的数据，列表的大小是 $18$ 行 $480$ 列，对应着 $18$ 个监测参数值和一个月的 $480$ 个小时。

#注意，由于range(n)函数生成的是 0 至 n-1的数，所以下面第 n 个月或第 n 天的索引是 n-1 
data_month={}   #创建一个空字典
for month in range(12):   #依次创建每个月的数据的列表 
    sample = np.empty([18,480])    #初始化一个列表
    for day in range(20):    #将每天的数据放到列表中对应的位置
        # 第 i 个月的数据是从第 i * 20 * 18 行开始的，第 j 天的数据是第 j * 18 行开始到 (j+1)*18 行结束。(PS: i 和 j 都是从 0 开始取，即第 1 个月对应 i = 0)
        # 第 j 天的数据应该存在 sample 的第 j * 24 列至 (j+1) * 24 列
        sample[:,day * 24 : (day+1) * 24] = raw_data[month * 18 * 20 + 18 * day : month * 18 * 20 + 18 * (day+1),:] 
    data_month[month] = sample   #将这个列表放入字典中，字典的 key 是月份，value 是对应的列表

Step3：重构数据得到 $\mathbf{X}$ 和 $\mathbf{y}$

  接下来是从上面的按月排好的数据中抽取出依次抽取出 9 个连续小时的数据作为 ${\mathbf{X}}_{train}$ 里的一条数据，后面的第 10 小时的 $\mathrm{P}\mathrm{M}2.5$ 取值作为对应的 ${\mathbf{y}}_{train}$ 的一个标签。例如对于 1 月，取 $0-8h$ 的 $18$ 种监测数据作为 ${\mathbf{X}}_{train}$ 的第一条数据 ${\mathbf{x}}_{train}^{(1)}$ ，第 $9h$ 的 $\mathrm{P}\mathrm{M}2.5$ 作为对应的标签 ${\mathbf{y}}_{train}^{(1)}$ 。第 $1-9h$ 的 $18$ 种监测数据作为 ${\mathbf{X}}_{train}$ 的第二条数据 ${\mathbf{x}}_{train}^{(2)}$ ，第 $10h$ 的 $\mathrm{P}\mathrm{M}2.5$ 作为对应的标签 ${\mathbf{y}}_{train}^{(2)}$ 。由于 $9$ 个小时的检测数据在上面的列表中是一个 ${\mathbb{R}}^{18\times 9}$ 的矩阵，所以还要将这个矩阵拉成一个 ${\mathbb{R}}^{1\times 162}$ 的行向量，可以调用 reshape(1,-1)来实现，拉取的结果就是原来第 $i$ 行第 $j$ 列的元素变到了新的行向量的第 $(i-1)\times 18+j$ 列上去了。

X = np.empty([12 * 471, 18 * 9], dtype = float)  #预分配 x 的空间，共有 12 个月乘 471 条数据，每条数据是 18(种) * 9(h) 维的
y = np.empty([12 * 471, 1], dtype = float)  #预分配 y 的空间，是 12*471 行 1 列的列向量
for month in range(12):  #共12个月
    for number in range(471):  #每个月有471条数据
        X[month * 471 + number, :] = data_month[month][:, number : number + 9].reshape(1,-1)
        y[month * 471 + number] = data_month[month][9, number + 9]  #pm2.5的数据在第 9 行	

  至此我们就完成了数据集的重构，得到了 $\mathbf{X}$ 和他的标签 $\mathbf{y}$ 。可以使用下面的代码将 $\mathbf{X}$ 和 $\mathbf{y}$ 保存成一个 $.\mathrm{c}\mathrm{s}\mathrm{v}$ 文件（下面这部分代码无需执行，只是贴出来我是如何保存的）

import csv 
with open('r_train_X.csv', mode='w', newline='') as r_train_X:   #保存 X
    csv_writer = csv.writer(r_train_X)
    for i in range(12 * 471):
        row = X[i][:]
        csv_writer.writerow(row)
with open('r_train_y.csv', mode='w', newline='') as r_train_y:   #保存 y
    csv_writer = csv.writer(r_train_y)
    for i in range(12 * 471):
        row = y[i]
        csv_writer.writerow(row)

直接读取重构好的数据

  如果觉得上面的从 $\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}.\mathrm{c}\mathrm{s}\mathrm{v}$ 文件重构数据得到 $\mathbf{X}$ 和 $\mathbf{y}$ 太麻烦了，操作不好。可以下载我已经重构好的数据文件 $\mathrm{r}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\_\mathrm{X}.\mathrm{c}\mathrm{s}\mathrm{v}$ 和 $\mathrm{r}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\_\mathrm{y}.\mathrm{c}\mathrm{s}\mathrm{v}$ 文件，再运行以下代码读取文件就可以得到上面所有步骤的结果。

import pandas as pd   
import numpy as np

data_X = pd.read_csv('./r_train_X.csv', encoding = 'big5')  #r_train_X.csv文件要与python代码文件(.py或.ipynb)在同一个目录下才行
data_y = pd.read_csv('./r_train_y.csv', encoding = 'big5')  #r_train_y.csv文件要与python代码文件(.py或.ipynb)在同一个目录下才行
X = data_X.to_numpy()  #将数据转换成numpy数据
y = data_y.to_numpy()  #将数据转换成numpy数据

线性回归模型代码实现（Python）

  下面正式开始实现一个线性回归模型。

Step1：数据预处理

  其实在上面说从 $\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}.\mathrm{c}\mathrm{s}\mathrm{v}$ 文件中重构出的 $\mathbf{X}$ 和 $\mathbf{y}$ 是训练集是不太准确的，其实还应该划分成训练集和验证集这两个部分，验证集用于调整模型的参数，所以在上面所说的所有的都应该只是 $\mathbf{X}$ ，下面划分后的才是真正用于训练的训练集。我采用的划分的比例为 $7:3$ ，划分后的数据分别记为 ${\mathbf{X}}_{train}、{\mathbf{y}}_{train}$ 和 ${\mathbf{X}}_{val}、{\mathbf{y}}_{val}$ ，划分的方式为随机选取。下面是划分的代码实现：

import math
X_rows, X_cols = X.shape   #获取 X 的行数 X_rows 和列数 X_cols
train_index = np.random.choice(np.arange(X_rows), size=math.floor(X_rows*0.7), replace=False)  #从 0 至 X_rows - 1 中随机挑选 (X_rows * 0.7) 向下取整个数，并且这些数不能重复，作为挑选出的训练集的索引
val_index = np.delete(np.arange(X_rows),train_index)  #剩下没被选到的作为验证集的索引
X_train = X[train_index,:]    #从 X 中取出被抽取到的训练集数据保存在 X_train 中
y_train = y[train_index,:]    #保存 X_train 对应的标签 y_train
X_val = X[val_index,:]        #从 X 中取出被抽取到的验证集数据保存在 X_val 中
y_val = y[val_index,:]        #保存 X_val 对应的标签 y_val

  如果留心观察数据集中各个属性的取值范围就可以发现取值范围差别是比较大的，因为后面寻找损失函数最小值的时候使用的是梯度下降，所以需要将数据缩放到相近的取值范围中，缩放的方式采用减去均值除以方差，即：
$$x_j^{(i)}\leftarrow \frac{x_j^{(i)}-{\mu }_j}{{\sigma }_j}$$
  并且缩放每个属性使用的均值 ${\mu }_j$ 和 方差 ${\sigma }_j$ 需要保存下来，因为在模型训练完在实际中进行预测时，也需对数据进行同样的缩放在输入到模型中。缩放的代码如下：

mean_X = np.mean(X_train, axis = 0)   #计算每个属性的均值
std_X = np.std(X_train, axis = 0)     #计算每个属性的方差
X_train = X_train - mean_X            #这里用到了numpy的传播特性来让每个数据减去对应的属性
X_train = X_train / std_X             #这里同样适用了numpy的传播特性来让每个数据除以对应的方差
X_val = X_val - mean_X                #在这预先对 X_val也做相同的缩放，这个也可以放在后面调模型参数的时候再做
X_val = X_val / std_X

Step2：模型训练

  下面使用梯度下降训练一个线性回归模型，首先我们在所有的数据的第一列的前面再加一列全为 $1$ 的属性记为 $x_0$ ，系数向量中也增加一个 ${\theta }_0$ ，这样就可以简化 $y=wx+b$ 成向量运算 $y=X\theta$ 。线性回归模型的假设函数 $h_{\theta }(x)$ 和损失函数 $J(\theta )$ 公式如下：
$$\begin{gathered} h_{\mathbf{\theta }}(\mathbf{x})={\mathbf{x}}^T\mathbf{\theta } \\ J(\mathbf{\theta })=\frac{1}{2m}\sum _{i=1}^m(h_{\mathbf{\theta }}({\mathbf{x}}^{(i)})-y^{(i)}{)}^2 \end{gathered}$$
  梯度的计算和系数的更新方式如下：
$$\frac{\partial J(\theta )}{\partial {\theta }_j}=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_j$$

$$\mathbf{\theta }:=\mathbf{\theta }-\eta \nabla J(\mathbf{\theta })$$

  下面是代码实现：

m, dim = X_train.shape   #获取 X 的行数 X_rows 和列数 X_cols
m_val = len(X_val)
dim = dim + 1  #维度加 1 \
X_train = np.concatenate((np.ones([len(X_train),1]),X_train), axis = 1).astype(float)  #在X_train前面加了一列1
X_val = np.concatenate((np.ones([len(X_val),1]),X_val), axis = 1).astype(float)  #在X_val前面加了一列1

theta = np.zeros([dim, 1 ])  #初始化系数向量
iteration = 1000     #预设迭代次数
learning_rate = 0.01    #预设学习率
loss_history = np.zeros([iteration, 1])    #用来记录每次梯度下降的损失值，用来画图观察梯度下降收敛的快慢
#loss_val_history = np.zeros([iteration, 1])   #用来记录在验证集上的损失值变化
for ite in range(iteration):
    h = np.dot(X_train, theta)   #计算每个样本输入到模型的输出
    h_y = h - y_train            #计算模型的预测和标签的差值，便于后面计算损失值和梯度
    loss = (1 / (2 * m)) * np.sum( np.power(h_y, 2) )   #计算损失值
    loss_history[ite] = loss
    
    # 记录在验证集上的损失值变化
    #pre = np.dot(X_val, theta)
    #loss_val_history[ite] = (1 / (2 * m_val)) * np.sum( np.power(pre - y_val, 2) )
    
    if (ite % 100 == 0):         #每 100 次迭代输出一次 loss 值
        print(str(ite) + " 's loss is: " + str(loss)) 
    gradient = (1 / m) * np.dot(X_train.T, h_y)    #计算梯度
    theta = theta - learning_rate * gradient       #更新梯度

  使用下面的代码画出损失值随迭代次数变化的曲线：

import matplotlib.pyplot as plt
itera = np.arange(iteration)
plt.plot(itera, loss_history)
#plt.plot(itera, loss_val_history)  #画出在验证集上损失值变化曲线
plt.show

  在上面的参数下画出来的曲线是下图这样的，可以看到基本收敛到了比较平稳的位置，应该是收敛到了最小值附近，将上面的部分代码的注释去掉，还以得到训练过程中在验证集上表现变化。

Step3：在验证集上测试

 在验证集上测试的过程和直接使用模型进行预测的过程是想通的，不过由于在前期已经对验证集的数据完成了正规化预处理，所以使用测试集数据继续回归预测时需要先进行前面的预处理步骤。由于 $\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}.\mathrm{c}\mathrm{s}\mathrm{v}$ 中没有标签结果，所以这里不在测试集上进行测试了，你也可以将数据集划分为训练集、验证集和测试集三个来自己生产一个测试集，不过用于训练的数据会变少了。

pre = np.dot(X_val, theta)   #计算预测结果
loss = (1 / (2 * m_val)) * np.sum( np.power(pre - y_val, 2) )   #计算损失值  

Step4：模型的改进

   在训练集上，样本的平均损失值大概是 16 左右，这个数还是非常大的，说明线性回归模型在训练集上的表现其实不是很好，需要使用更复杂的模型，比如不仅使用一次函数，而使用多项式函数。下面尝试使用 $2$ 次多项式函数。为了减少参数，不考虑 $x_i{x}_j$ 这种情况，仅考虑 $x_i^2$ 这种。

下载链接

① 原始数据集 $\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}.\mathrm{c}\mathrm{s}\mathrm{v}$ 和 $\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}.\mathrm{c}\mathrm{s}\mathrm{v}$ 文件

链接：https://pan.baidu.com/s/1_OtBbNJS2v3vIPbiaulNgQ
提取码：8i0u

② 重构后的数据集 $\mathrm{r}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\_\mathrm{X}.\mathrm{c}\mathrm{s}\mathrm{v}$ 和 $\mathrm{r}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\_\mathrm{y}.\mathrm{c}\mathrm{s}\mathrm{v}$

链接：https://pan.baidu.com/s/19A0YrxyBbmeRw9tAHDbaqA
提取码：6o5p

③ Matlab 的 $.\mathrm{m}\mathrm{a}\mathrm{t}$ 文件如有需要的人我再补