前言
李宏毅2020机器学习的第一个作业是做PM2.5预测，模型用线性回归（linear regression），优化方法用梯度下降（gradient descent），要求手写。老师给了训练集（training set）和测试集（testing set）和参考代码，本文是在Pycharm上的实现的，pycharm上的环境搭建见另一篇博客，这个课程后期也用Pytorch。需要完整pycharm工程代码、数据集和测试结果等资料点击此处，没有积分的，可以私聊我，不过可能需要等我看到。理论见这篇博客
Pytorch版见：深度之眼pytorch打卡（四）| 台大李宏毅机器学习 2020作业（一）：线性回归，实现多因素作用下的PM2.5预测（Pytorch版）
原理分析见：台大李宏毅 机器学习 2020学习笔记（二）：回归与过拟合
数据集
- 训练集
在网上下载的数据集编码方式是“big5”，繁体字的编码方式，下载后直接用excell打开会出现乱码。可以在Pycharm中打开，不过要先选择编码方式为“Big5-HKSCS”，再打开文件才不会出现乱码。正常打开文件之后，再选择编码方式为GBK，点击convert，则原来的文件就会变成GBK编码方式。
训练集包含2014年12个月中每个月取20天，每天24小时的环境监测数据，每个小时的数据都有18各特征，包括AMB_TEMP（温度）, CH4, CO, NHMC, NO, NO2, NOx, O3, PM10, PM2.5, RAINFALL（降雨量）, RH, SO2, THC, WD_HR, WIND_DIREC, WIND_SPEED, WS_HR。部分训练集数据如下图。
- 测试集
测试集有训练集类似的格式，共240个id，每个id有9小时的数据，每组数据有18个特征。要求根据每个id前9个小时的所有数据来预测第10个小时的PM2.5值。
数据处理
主要参照老师给的参考文档
- 加载训练集并预处理
新建一个Pycharm工程，并在工程中建一个data文件夹用于存放数据。训练集中前三列为非数据内容，需要去掉，并且RAINFALL特征中的NR，即Not Rain 需要用数字0替换。
iloc：根据行列号来取数据。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
data = pd.read_csv('./data/train.csv', encoding='big5')
data = data.iloc[:, 3:] # 去掉前三列，即非数据列
data[data == 'NR'] = 0 # 将rainfall里的NR换成0，即不下雨为0
raw_data = data.to_numpy
# print(raw_data)- 训练集分组
为了便于处理，先将每个月的资料由360*24的数据，转换成18*480的数据，如下图所示（图引自老师给的参考文档）。代码先将每个月的数据放到一个18*480的二维数组中，然后再将二维数组整体放入一个字典。最后生成的字典，以月份为索引，然后每个索引对应一个18*480的二维数组。
month_data = {}
for month in range(12):
sample = np.empty([18, 480])
for day in range(20):
sample[:, day*24:(day+1)*24] = raw_data[18*(day+20*month):18*(day+1+20*month), :] # 从一天的开始取到一天的结尾
month_data[month] = sample
# print(month_data)由于测试要求是根据前九个小时的数据推测出第十个小时的PM2.5值，所以再把10个小时的数据分为一组，其中前9个小时的18*9个特征作为输入x，第10个小时的第10个特征Pm2.5值作为输出y。采用步长为1的移动取值方式依次取10个值，如下图所示（图引自老师给的参考文档）则一月可取出（480-10）/1 +1=471组数据。那么一年就有12*471组数据。
reshape(1, -1)：转换成一个行向量。
reshape(-1, 1) ：转换成一个列向量。
参考文档在进行移位取值时用天和时间，这使的取值的参数变得很复杂，既然已经将一个月的数据合并，那么天与小时都已经是次要的信息。因此，此处对于1个月，只需要一个index，从0到471取到就可以了。
x = np.empty([12*471, 18*9], dtype=float)
y = np.empty([12*471, 1], dtype=float)
for month in range(12):
# 注释掉的部分为老师给的参考代码
# for day in range(20):
# for hour in range(24):
# if hour > 14 and day == 19: # 在每个月的结尾处，hour大于14后就没有一组数据了，即大于471小时
# continue
# x[month*471 + day*24 + hour, :] = month_data[month][:, day*20 + hour:day*20+hour+9].reshape(1, -1)
# # 以month值为索引，将每一组前9个小时的9*18的二维数组转换成一个9*18维的行向量，并赋值给x中的行
# y[month*471 + day*24 + hour, :] = month_data[month][9, day*20+hour+9]
# # 以month值为索引，将每一组第10个小时第10个特征，即PM2.5值给y的行
for index in range(471):
x[index+471*month, :] = month_data[month][:, index:index+9].reshape(1, -1)
y[index + 471 * month, :] = month_data[month][9, index + 9]
print(x[10])
print(y[10])对照x与训练集2014年1月1日10-18点内容，和y与19点的PM2.5值，可知是吻合的，代码没有问题。
-
训练集标准化（Normalize）
数据标准化之后模型更容易收敛，也更快收敛，方法就是让其减去均值再除以方差。
np.mean中，axis=0表示列均值，axis=1表示行均值。
mean_x = np.mean(x, axis=0) # 列均值
std_x = np.std(x, axis=0)
for i in range(len(x)):
for j in range(len(x[0])):
if std_x[j]!=0:
x[i][j] = (x[i][j]-mean_x[j])/std_x[j] # 每个元素减去其所在列的均值，再除以标准差
模型训练
模型：线性函数
输入：一次18*9个数据
参数：18*9个权重+一个偏置
损失函数：均方根误差（Root Mean Square Error）
优化算法：梯度下降、Adagrad
np.concatenate((a,b),axis=1) ： 进行行拼接，拼接结果是[a,b]
np.dot(a,b)：a与b点乘
np.power(a,b)：a的b次方根
- 梯度下降法（BGD）优化
dim = 18*9 + 1
w = np.zeros([dim, 1])
x = np.concatenate((np.ones([12*471, 1]), x), axis=1).astype(float) # 为1那一列是bias前的系数
N = 12*471
learning_rate = 0.03
iter_time = 1000
loss = np.zeros([iter_time, 1])
for i in range(iter_time):
y1 = np.dot(x, w)
loss[i] = np.sqrt(np.sum(np.power(y1-y, 2))/N)
gradient = 2*np.dot(x.transpose(), y1-y)/N # 比原文，笔者在此处多除了个N，这样可以用更大的学习率
w = w - learning_rate*gradient
plt.plot(loss)
plt.show() 
从左到右依次是0.03,0.05和0.001的学习率，迭代了1000次，较高的学习率无法收敛，较低的学习率收敛慢。
- Adagrad优化
比起BGD，这种优化方法，在更新参数时，除了乘学习率外，还除了过往所有梯度平方和的平方根，可动态的改变学习率。并且不会因为某时梯度太大而导致模型发散。
dim = 18*9 + 1
w = np.zeros([dim, 1])
x = np.concatenate((np.ones([12*471, 1]), x), axis=1).astype(float) # 为1那一列是bias前的系数
N = 12*471
learning_rate = 0.03
iter_time = 1000
loss = np.zeros([iter_time, 1])
adagrad = np.zeros([dim, 1])
eps = 0.000000001
for i in range(iter_time):
y1 = np.dot(x, w)
loss[i] = np.sqrt(np.sum(np.power(y1-y, 2))/N)
gradient = 2*np.dot(x.transpose(), y1-y)/N # 比原文，笔者在此处多除了个N，这样可以用更大的学习率
adagrad += gradient**2
w = w - learning_rate*gradient/np.sqrt(adagrad+eps) # eps为了防止分母为零的一个很小的数
plt.plot(loss)
plt.show()
np.save('weight.npy', w) 
从左到右学习率是0.03,1，和100。
预测
由老师给的测试集，是没有第十个小时的，即相当于是直接预测，而不是测试。如果要做一下预测，可以从训练集从拿一个月的数据出来做测试集。
ValueError: could not broadcast input array from shape (153) into shape (162)
在处理测试集的时候出现这个错，是因为测试集没有列号，第一行就是数据，导致pd.read_csv将第一行数据当成了列号，因此数据少了一行。到最后就出现了维度对不上的错误。如图所示，本应该是4320*9。
test_data = pd.read_csv('./data/test.csv', encoding='big5', header=None) #加个header=None即可解决。 
- 预测
数据的处理方式都与训练集相同，包括点乘生成y。由于预测的结果有负值，在输出时我将负值令成了0。
import csv
import pandas as pd
import numpy as np
test_data = pd.read_csv('./data/test.csv', encoding='big5', header=None)
test_data = test_data.iloc[:, 2:]
test_data[test_data == 'NR'] = 0
test_raw_data = test_data.to_numpy()
test_x = np.empty([240, 18*9], dtype=float)
# print(test_raw_data)
for i in range(240):
test_x[i, :] = test_raw_data[18 * i: 18 * (i + 1), :].reshape(1, -1)
# print(test_x[0])
test_mean = np.mean(test_x, axis=0)
test_std = np.std(test_x, axis=0)
for i in range(len(test_x)):
for j in range(len(test_x[0])):
if test_std[j] != 0:
test_x[i][j] = (test_x[i][j] - test_mean[j])/test_std[j]
# print(test_x[0])
test_x = np.concatenate((np.ones([240, 1]), test_x),axis=1).astype(float)
test_y = np.empty([240, 1], dtype=float)
w = np.load('weight.npy')
test_y = np.dot(test_x, w)
for i in range(240):
if int(test_y[i][0]) > 0:
print('id:', i, ' PM2.5:', int(test_y[i][0]))
else:
print('id:', i, ' PM2.5:', 0)
- 输出到CSV文件
with open('predict.csv', mode = 'w', newline='') as file:
csv_writer = csv.writer(file)
header = ['id', 'PM2.5']
csv_writer.writerow(header)
for i in range(240):
if int(test_y[i][0]) > 0:
row = ['id_'+str(i), str(int(test_y[i][0]))]
else:
row = ['id_'+str(i), '0']
csv_writer.writerow(row)
最低0.47元/天 解锁文章
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
