线性回归预测PM2.5----台大李宏毅机器学习做业1(HW1)

1、做业说明

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

训练集介绍：

　　(1)、CSV文件，包含台湾丰原地区240天的气象观测资料(取每月前20天的数据作训练集，12月X20天=240天，每个月后10天数据用于测试，对学生不可见);

　　(2)、天天的监测时间点为0时，1时......到23时，共24个时间节点;

　　(3)、天天的检测指标包括CO、NO、PM2.五、PM10等气体浓度，是否降雨、刮风等气象信息，共计18项；

　　(4)、数据集地址：https://pan.baidu.com/s/1o2Yx42dZBJZFZqCa5y3WzQ，提取码：qgtm。 

 用excel打开，繁体字会出现乱码：

用notepad++能够正常打开：

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替换指定元素，将空数据填充为0
    df = df.replace(['NR'], [0.0])
    # astype() 转换array中元素数据类型
    array = np.array(df).astype(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.5
            x_list.append(mat)
            y_list.append(label)
    x = np.array(x_list)
    y = np.array(y_list)
    
    return x, y, array

2.2 模型创建

　　若是对相关领域比较熟悉的话，能够根据PM2.5与PM十、SO、NO的浓度关系选择合适的模型。

　　若是对数据比较敏感的话，能够从数据中发现规律并以此为依据创建模型。

　　不过笔者对气象领域并不熟悉，对数据也不够敏感，只能采用最简单、最low的线性回归模型。不过既然是做业嘛，就应该容许学生随意发挥，不见得就存在标准答案。

2.2.1 回归模型

　　采用最普通的线性回归模型，并无用上训练集中全部的数据，只用到了每一个数据帧样本中的9个PM2.5含量值：

为对应数据帧中第i个PM2.5含量，为其对应的权重值，为偏置项，为该数据帧样本的预测结果。                                   

2.2.2 损失函数

　　用预测值与label之间的平均欧式距离来衡量预测的准确程度，并充当损失函数（这里的损失指的是平均损失；乘1/2是为了在后续求梯度过程当中保证梯度项系数为1，方便计算）：

 

为第n个label，为第n个数据帧的预测结果，为参加训练的数据帧样本个数。

　　为了防止过拟合，加入正则项：

 

为正则项，为正则项系数。

2.2.3 梯度更新

　　梯度计算：需明确此时的目标是使Loss最小，而可优化的参数为权重w和偏置值b，所以须要求Loss在w上的偏微分和Loss在b上的偏微分。

　　计算出梯度后，经过梯度降低法实现参数更新。

为权重w更新时的学习率，为偏置b更新时的学习率。

2.2.3 学习率更新

　　为了在不影响模型效果的前提下提升学习速度，能够对学习率进行实时更新：即让学习率的值在学习初期较大，以后逐渐减少。这里采用比较经典的adagrad算法来更新学习率。

为更新后的学习率，为更新前的学习率。为在此以前全部梯度平方和的二次根。

# 更新参数，训练模型
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 = 0
        w_g = np.zeros(9)
        # 在全部数据上计算Loss_label的梯度
        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]) 
        # 求平均    
        b_g /= 3200
        w_g /= 3200
        #  加上Loss_regularization在w上的梯度
        for m in range(9):
            w_g[m] += reg_rate * weights[m]
        
        # adagrad
        bg2_sum += b_g**2
        wg2_sum += w_g**2
        # 更新权重和偏置
        bias -= learning_rate/bg2_sum**0.5 * b_g
        weights -= learning_rate/wg2_sum**0.5 * w_g

    return weights, bias

 3、代码分享与结果分析

3.1 源代码

 1 import pandas as pd
 2 import numpy as np
 3 
 4 # 数据预处理
 5 def dataProcess(df):
 6     x_list, y_list = [], []
 7     # df替换指定元素，将空数据填充为0
 8     df = df.replace(['NR'], [0.0])
 9     # astype() 转换array中元素数据类型
10     array = np.array(df).astype(float)
11     # 将数据集拆分为多个数据帧
12     for i in range(0, 4320, 18):
13         for j in range(24-9):
14             mat = array[i:i+18, j:j+9]
15             label = array[i+9, j+9] # 第10行是PM2.5
16             x_list.append(mat)
17             y_list.append(label)
18     x = np.array(x_list)
19     y = np.array(y_list)
20 
21     '''
22     # 将每行数据都scale到0到1的范围内，有利于梯度降低，但经尝试发现效果并很差
23     for i in range(18):
24         if(np.max(x[:, i, :]) != 0):
25             x[: , i, :] /= np.max(x[:, i, :])
26     '''
27     return x, y, array
28 
29 # 更新参数，训练模型
30 def train(x_train, y_train, epoch):
31     bias = 0 # 偏置值初始化
32     weights = np.ones(9) # 权重初始化
33     learning_rate = 1 # 初始学习率
34     reg_rate = 0.001 # 正则项系数
35     bg2_sum = 0 # 用于存放偏置值的梯度平方和
36     wg2_sum = np.zeros(9) # 用于存放权重的梯度平方和
37 
38     for i in range(epoch):
39         b_g = 0
40         w_g = np.zeros(9)
41         # 在全部数据上计算Loss_label的梯度
42         for j in range(3200):
43             b_g += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-1)
44             for k in range(9):
45                 w_g[k] += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-x_train[j, 9, k]) 
46         # 求平均    
47         b_g /= 3200
48         w_g /= 3200
49         #  加上Loss_regularization在w上的梯度
50         for m in range(9):
51             w_g[m] += reg_rate * weights[m]
52         
53         # adagrad
54         bg2_sum += b_g**2
55         wg2_sum += w_g**2
56         # 更新权重和偏置
57         bias -= learning_rate/bg2_sum**0.5 * b_g
58         weights -= learning_rate/wg2_sum**0.5 * w_g
59 
60         # 每训练200轮，输出一次在训练集上的损失
61         if i%200 == 0:
62             loss = 0
63             for j in range(3200):
64                 loss += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias)**2
65             print('after {} epochs, the loss on train data is:'.format(i), loss/3200)
66 
67     return weights, bias
68 
69 # 验证模型效果
70 def validate(x_val, y_val, weights, bias):
71     loss = 0
72     for i in range(400):
73         loss += (y_val[i] - weights.dot(x_val[i, 9, :]) - bias)**2
74     return loss / 400
75 
76 def main():
77     # 从csv中读取有用的信息
78     # 因为你们获取数据集的渠道不一样，因此数据集的编码格式可能不一样
79     # 若读取失败，可在参数栏中加入encoding = 'gb18030'
80     df = pd.read_csv('train.csv', usecols=range(3,27))
81     x, y, _ = dataProcess(df)
82     #划分训练集与验证集
83     x_train, y_train = x[0:3200], y[0:3200]
84     x_val, y_val = x[3200:3600], y[3200:3600]
85     epoch = 2000 # 训练轮数
86     # 开始训练
87     w, b = train(x_train, y_train, epoch)
88     # 在验证集上看效果
89     loss = validate(x_val, y_val, w, b)
90     print('The loss on val data is:', loss)
91 
92 if __name__ == '__main__':
93     main()

View Code

3.1 结果展现

　　能够看出，模型在验证集上的损失为40左右，即预测值与label之间的平均差别在6到7之间，因而可知，模型的总体效果仍是比较差的。

3.3 模型改进的方向

　　(1)在从csv文件中提取数据帧和label时，本文以天为单位，天天分割出15个数据帧和15个label。事实上，时间是连续的，能够将每个月的20天首尾链接，再从其中分割数据帧和label，可以使数据帧样本数量大大提高，可能会使模型效果更优。

　　(2)在构建模型时，应充分考虑PM2.5与其余大气成分之间的关系，构建更合理的模型。

　　(3)分割训练集和验证集时，应该按照比例随机抽取数据帧做为训练集和验证集，而不是像本文那样简单地把前3200个数据样本做为训练集，后400个做为验证集。

 

参考资料：

　　李宏毅老师机器学习课程视频：https://www.bilibili.com/video/av10590361

　　李宏毅老师机器学习课程讲义资料：http://speech.ee.ntu.edu.tw/~tlkagk/courses_ML17_2.html

　　邱锡鹏老师《神经网络与深度学习》： https://nndl.github.io/