一, 说明
给定训练集 train.CSV, 要求根据前 9 个小时的空气监测情况预测第 10 个小时的 PM2.5 含量.
训练集介绍:
(1),CSV 文件, 包含台湾丰原地区 240 天的气象观测资料 (取每个月前 20 天的数据做训练集, 12 月 X20 天 = 240 天, 每月后 10 天数据用于测试, 对学生不可见);
(2), 每天的监测时间点为 0 时, 1 时...... 到 23 时, 共 24 个时间节点;
(3), 每天的检测指标包括 CO,NO,PM2.5,PM10 等气体浓度, 是否降雨, 刮风等气象信息, 共计 18 项;
二, 数据处理
根据要求, 要用前 9 个小时的数据, 来预测第 10 个小时的 PM2.5 值. 一笔训练数据如下图所示:
数据中存在一定量的空数据 NR, 且多存在于 RAINFALL 一项. 对于空数据, 常规的处理方法无非就是删除法和补全法两种.
RAINFALL 表示当天对应时间点是否降雨, 因此可以采用补全法处理空数据: 将空数据 NR 全部补为 0 即可.
- # 将 NR 替换成 0
- data = data.replace(['NR'], [0.0])
我们先将数据进行预处理, 得到 每笔训练数据 和对应的 结果 label.
1. 由于每个月前 20 天的数据检测是连续 24 小时进行的, 所以为了得到多笔数据, 先将每个月 20 天数据 "连起来", 如下图示:
每个月的数据就是 18 行 480(24*20) 列, 一共 12 个月的数据.
- # 将每个月 20 天数据连成一大行
- month_data = []
- for month in range(12):
- # 每个月的数据
- sub_data = np.empty([18, 20*24])
- for day in range(20):
- # 每一天的数据
- sub_data[:, day*24:(day+1)*24] = data[(month*18*20+day*18):(month*18*20+(day+1)*18), :]
- month_data.append(sub_data)
2. 对于连续的 10 个小时, 可以取一笔 9 小时训练数据 (输入) 和 第 10 小时对应 PM2.5 值 (结果).
每个月 20 天, 20*24=480 小时, 480-9=471, 每个月可以取 471 笔数据.
- # 将每个月中 20 天, 相邻 9 个小时生成一笔数据, 第 10 个小时的 pm2.5 值, 生成一个 label
- for i in range(12):
- sub_data = month_data[i]
- for j in range(20*24-9):
- # 相邻 9 小时的数据
- x_list.append(sub_data[:, j:j+9])
- # 第 10 小时的 pm2.5
- y_list.append(sub_data[9, j+9])
完整数据处理代码:
- def data_process(data):
- x_list , y_list = [], []
- # 将 NR 替换成 0
- data = data.replace(['NR'], [0.0])
- # astype() 转换为 float
- data = np.array(data).astype(float)
- # 将每个月 20 天数据连成一大行
- month_data = []
- for month in range(12):
- # 每个月的数据
- sub_data = np.empty([18, 20*24])
- for day in range(20):
- # 每一天的数据
- sub_data[:, day*24:(day+1)*24] = data[(month*18*20+day*18):(month*18*20+(day+1)*18), :]
- month_data.append(sub_data)
- # 将每个月中 20 天, 相邻 9 个小时生成一笔数据, 第 10 个小时的 pm2.5 值, 生成一个 label
- for i in range(12):
- sub_data = month_data[i]
- for j in range(20*24-9):
- # 相邻 9 小时的数据
- x_list.append(sub_data[:, j:j+9])
- # 第 10 小时的 pm2.5
- y_list.append(sub_data[9, j+9])
- x = np.array(x_list)
- y = np.array(y_list)
- return x, y, month_data
- View Code
二, 模型建立
如题所说, 使用的是最简单的线性回归模型, 作为课程作业没有那么难, 但也从中学到不少东西.
2.1 线性回归模型
如果把 b 作为 w0, 加到权值向量前面, 可以得到向量运算的形式, 如下:
h(X) = WTX #W 为权值, X 为输入.
2.2 损失函数
用预测值与 label 之间的平均欧式距离来衡量预测的准确程度, 并充当损失函数.
这里的损失指的是平均损失; 乘 1/2 是为了在后续求梯度过程中保证梯度项系数为 1, 方便计算.
为了防止过拟合, 加入正则项:
完整的损失函数:
2.3 梯度下降
对参数 w 和 b 求偏导:
对参数进行更新:
ηw,ηb 为学习率.
2.4 学习率更新
为了在不影响模型效果的前提下提高学习速度, 可以对学习率进行实时更新: 即让学习率的值在学习初期较大, 之后逐渐减小.
这里采用比较经典的 adagrad 算法来更新学习率:
根号下为梯度的累加值.
2.5 矩阵加速计算推导
因为 python 中使用矩阵的计算速度非常快, 远远快于循环计算, 所以这里我们推导一下利用矩阵计算梯度值的写法.
具体可以参照这篇博客:
前文提到过, 如果把 b 作为 w0, 加到权值向量前面可以得到如下形式:
h(X) = WTX # W 为权值, X 为输入. W = [b, w0, w1, ...]
对于损失函数的转化改造:
对于
1. 将平方求和改造成向量模的平方:
假设我们有个向量
所以有
将平方项展开:
对 W 进行求偏导:
与上面展开对应
所以我们要求的梯度就是:
* 这里将分子的 2 去掉是因为与的损失函数分母抵消.
- # 计算梯度 W = X 转置.(XW-Y)
- w_1 = np.dot(X.transpose(), X.dot(W)-y_train)
三, 训练模型
3.0 数据转化
将训练数据分成两部分 (8:2), 一部分用来训练, 一部分用来验证效果.
- # 8:2 cross validation
- x_train = x[:(int)(x.shape[0]*0.8)]
- y_train = y[:(int)(x.shape[0]*0.8)]
- x_val = x[(int)(x.shape[0]*0.8+0.5):]
- y_val = y[(int)(y.shape[0]*0.8+0.5):]
由于参数太多, 也可以取其中的几类进行训练, 比如下文中将选取 NO,NO2,NOx,O3,PM10,PM2.5 作为输入.
其中的一笔数据如下:
首先将每笔数据的输入转化成一行, 并在前面加上 1, 对应于 bias 项.
- # 定义参数 b,w b 作为 w0
- W = np.ones(1+9*6)
- # 将训练数据转化成 每一笔数据一行, 并且前面添加 1, 作为 b 的权值 [[1, ...], [1, ...],...,[1, ...]]
- X = np.empty([n, W.size-1])
- for i in range(n):
- X[i] = x_train[i][4:10].reshape(1, -1)
- # 添加 1
- X = np.concatenate((np.ones([n, 1]), X), axis=1)
3.1 训练函数
完整的训练函数代码如下, 具体请看注释:
- def train(x_train, y_train, times):
- # 定义参数 b,w b 作为 w0
- W = np.ones(1+9*6)
- # 多少笔数据
- n = y_train.size
- # 学习率
- learning_rate = 100
- # 正则项大小
- reg_rate = 0.011
- # 将训练数据转化成 每一笔数据一行, 并且前面添加 1, 作为 b 的权值 [[1, ...], [1, ...],...,[1, ...]]
- X = np.empty([n, W.size-1])
- for i in range(n):
- X[i] = x_train[i][4:10].reshape(1, -1)
- # 添加 1
- X = np.concatenate((np.ones([n, 1]), X), axis=1)
- # 累加正则项
- adagrad=0
- # 正则项的选择矩阵, 去掉 bias 部分
- reg_mat=np.concatenate((np.array([0]), np.ones([9*6,])), axis=0)
- for t in range(times):
- # 计算梯度 W = X 转置.(XW-Y)
- w_1 = np.dot(X.transpose(), X.dot(W)-y_train)
- # 加正则项
- w_1 += reg_rate * W * reg_mat
- # 正则项参数更新
- adagrad += sum(w_1**2)**0.5
- # 梯度下降
- W -= learning_rate/adagrad * w_1
- # 每 200 次迭代输出一次
- if t%200==0:
- loss = 0
- for j in range(n):
- loss += (y_train[j]-X[j].dot(W))**2
- print(t)
- print('times', loss/n)
- return W
3.2 验证
- def validate(x_val, y_val, w):
- n=y_val.size
- # 转化成一行, 并加一列 1
- X = np.empty([n, w.size - 1])
- for i in range(n):
- X[i] = x_val[i][4:10].reshape(1, -1)
- X = np.concatenate((np.ones([n, 1]), X), axis=1)
- loss=0
- # 计算 loss
- for j in range(n):
- loss += (y_val[j] - X[j].dot(W)) ** 2
- return loss/n
四, 结果分析
运行输出的结果看, loss 还是挺大的, 还有改进的空间.
改进思路:
1. 分割训练集和验证集时, 应该按照比例随机抽取数据帧作为训练集和验证集, 选取 loss 最小的模型.
2. 充分考虑其他参数对空气 PM2.5 的影响, 加入更加复杂的高次项.
五, 预测结果
对 test 集的数据进行结果预测
- ## 计算预测值 ##
- Y = X_test.dot(W)
- # 预测值写入
- data_test = np.array(data_test)
- data_test = np.concatenate((data_test, np.zeros([n, 1])), axis=1)
- for j in range(0, n, 18):
- data_test[j+9, 11] = int(Y[int(j/18)]+0.5)
为了方便查看, 将数据写回源文件格式.
- # 保存结果
- data_test = pd.DataFrame(data_test)
- data_test.to_csv('test_res.csv')
第一笔数据的预测值:
六, 程序代码
** 在项目根目录存放'train.csv','test.csv'
** 每次训练后会保存参数, 下次训练时请事先删除根目录文件'weight_2.npy'
- import pandas as pd
- import numpy as np
- def data_process(data):
- x_list , y_list = [], []
- # 将 NR 替换成 0
- data = data.replace(['NR'], [0.0])
- # astype() 转换为 float
- data = np.array(data).astype(float)
- # 将每个月 20 天数据连成一大行
- month_data = []
- for month in range(12):
- # 每个月的数据
- sub_data = np.empty([18, 20*24])
- for day in range(20):
- # 每一天的数据
- sub_data[:, day*24:(day+1)*24] = data[(month*18*20+day*18):(month*18*20+(day+1)*18), :]
- month_data.append(sub_data)
- # 将每个月中 20 天, 相邻 9 个小时生成一笔数据, 第 10 个小时的 pm2.5 值, 生成一个 label
- for i in range(12):
- sub_data = month_data[i]
- for j in range(20*24-9):
- # 相邻 9 小时的数据
- x_list.append(sub_data[:, j:j+9])
- # 第 10 小时的 pm2.5
- y_list.append(sub_data[9, j+9])
- x = np.array(x_list)
- y = np.array(y_list)
- return x, y, month_data
- def train(x_train, y_train, times):
- # 定义参数 b,w b 作为 w0
- W = np.ones(1+9*6)
- # 多少笔数据
- n = y_train.size
- # 学习率
- learning_rate = 100
- # 正则项大小
- reg_rate = 0.011
- # 将训练数据转化成 每一笔数据一行, 并且前面添加 1, 作为 b 的权值 [[1, ...], [1, ...],...,[1, ...]]
- X = np.empty([n, W.size-1])
- for i in range(n):
- X[i] = x_train[i][4:10].reshape(1, -1)
- # 添加 1
- X = np.concatenate((np.ones([n, 1]), X), axis=1)
- # data_X = pd.DataFrame(X)
- # data_X.to_csv('data.csv')
- adagrad=0
- # 正则项的选择矩阵, 去掉 bias 部分
- reg_mat=np.concatenate((np.array([0]), np.ones([9*6,])), axis=0)
- for t in range(times):
- # 计算梯度 W = X 转置.(XW-Y)
- w_1 = np.dot(X.transpose(), X.dot(W)-y_train)
- # 加正则项
- w_1 += reg_rate * W * reg_mat
- # 正则项参数更新
- adagrad += sum(w_1**2)**0.5
- # 梯度下降
- W -= learning_rate/adagrad * w_1
- # 每 200 次迭代输出一次
- if t%200==0:
- loss = 0
- for j in range(n):
- loss += (y_train[j]-X[j].dot(W))**2
- print('After', t,'times loss=', loss/n)
- return W
- def validate(x_val, y_val, w):
- n=y_val.size
- # 转化成一行, 并加一列 1
- X = np.empty([n, w.size - 1])
- for i in range(n):
- X[i] = x_val[i][4:10].reshape(1, -1)
- X = np.concatenate((np.ones([n, 1]), X), axis=1)
- loss=0
- # 计算 loss
- for j in range(n):
- loss += (y_val[j] - X[j].dot(W)) ** 2
- return loss/n
- if __name__ == '__main__':
- data = pd.read_csv('./train.csv', encoding='big5')
- # 去掉前三列
- data = data.iloc[:, 3:]
- [x, y, month_data] = data_process(data)
- # 8:2 cross validation
- x_train = x[:(int)(x.shape[0]*0.8)]
- y_train = y[:(int)(x.shape[0]*0.8)]
- x_val = x[(int)(x.shape[0]*0.8+0.5):]
- y_val = y[(int)(y.shape[0]*0.8+0.5):]
- try:
- W = np.load('weight_2.npy')
- except:
- # 迭代次数
- times = 10000
- W = train(x_train, y_train, times)
- np.save('weight_2.npy', W)
- ## 计算在 val 上的 loss ##
- loss = validate(x_val, y_val, W)
- print('validate loss=', loss)
- ## 在 test 上进行验证 ##
- # header=None 无表头读入
- data_test = pd.read_csv('./test.csv', header=None, encoding='big5')
- # 去掉前两列
- test = data_test.iloc[:, 2:]
- test = test.replace(['NR'], [0.0])
- #处理数据
- test = np.array(test).astype(float)
- [n, m] = test.shape
- # 读出参数值
- X_test = np.empty([int(n/18), 9*6])
- for i in range(0, n, 18):
- X_test[int(i/18), :] = test[i+4:i+10, :].reshape(1, -1)
- [n_test, m_test] = X_test.shape
- # 加一列 1
- X_test = np.concatenate((np.ones([n_test, 1]), X_test), axis=1)
- ## 计算预测值 ##
- Y = X_test.dot(W)
- # 预测值写入
- data_test = np.array(data_test)
- data_test = np.concatenate((data_test, np.zeros([n, 1])), axis=1)
- for j in range(0, n, 18):
- data_test[j+9, 11] = int(Y[int(j/18)]+0.5)
- # 保存结果
- data_test = pd.DataFrame(data_test)
- data_test.to_csv('test_res.csv')
- View Code
来源: https://www.cnblogs.com/Dawn-bin/p/12813991.html