最近开始学习机器学习，听同学介绍，选择了台大李宏毅老师的视频，https://www.bilibili.com/video/BV1JE411g7XF
老师布置了很多作业，这里记录下自己写作业的过程和一些参考资料
本篇是作业一，根据丰源站历史的空气数据，预测下一小时的pm2.5的值。
本次作业说明视频参见 https://www.bilibili.com/video/BV1gE411F7td
import numpy as np
import pandas as pd
本次作业要求只用numpy,pandas手写线性回归，所以只导入pandas和numpy两个包，
学习这两个包推荐一个gitbook， 利用python进行数据分析
另外这本书有本实体版的，是另外一个人翻译的，翻译的非常差，非常不建议购买，看这个gitbook版的就好。
本次作业代码参考 https://zhuanlan.zhihu.com/p/115335196
上面知乎专栏的大佬写的比较简单，我添加了更多的细节，便于新人理解。
data = pd.read_csv('data/train.csv', encoding='big5')
首先使用pandas读取训练文件，并做一些预处理，李宏毅老师是湾湾人，
用的是繁体，所以要指定下编码为big5，否则里面的中文为乱码。
看一下文件的结构，有27列，第一列为日期，第三列为检测项，后面的是0-23点每一点的数值。
然后每天有18行，即18个检测项，其中第10行是pm2.5。
一共12月*20天*18项=4320行
data.head(20)
然后前三项是没用的，我们用切片去掉。shape就变为了4320*24
数据中还有很多是NR，我们替换为0。
接下来就是numpy的处理范围了，使用to_numpy转换为矩阵。
data = data.iloc[:,3:]
data.replace('NR',0,inplace=True)
row_data = data.to_numpy()
row_data
参考上图，我们需要把数据shape做一个转换，把240天内相同的测量项放在同一行。
也就是从(12月*20天*18项),24小时变为(18项,(24小时*20天),12月
下面的两个循环就是在做转换，每个月的数据写在以月份为键名的字典里的。
month_data={}
for month in range(12):
sample = np.empty([18, 480])
for day in range(20):
sample[:,day*24:(day+1)*24] = row_data[18*(month*20+day):18*(month*20+(day+1)),:]
month_data[month] = sample
然后我们以10小时为一个data，把前9个小时的data作训练，第10个小时的值为target。
即0-9点的data来预测10点的pm2.5，第1-10点的和data来预测11点的pm2.5
最后一笔是471-479来预测480的pm2.5，所以一共有471笔data
shape从变为
y就是第10个小时的pm2.5的值，shape为
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 day == 19 and hour > 14:
continue
x[month * 471 + day * 24 + hour, :] = \
month_data[month][:,day * 24 + hour : day * 24 + hour + 9].reshape(1, -1)
y[month * 471 + day * 24 + hour, 0] = month_data[month][9, day * 24 + hour + 9]
x
接下来做Normalize，教学视频里说明了为什么要做Normalize，不知道的请复习视频
按照上面的公式，需要求平均值和标准差，然后x处理对应位置的值
mean_x = np.mean(x, axis = 0)
std_x = np.std(x, axis = 0) #18 * 9
for i in range(len(x)): #12 * 471
for j in range(len(x[0])): #18 * 9
if std_x[j] != 0:
x[i][j] = (x[i][j] - mean_x[j]) / std_x[j]
x
接下来就是最重要的部分了，我们开始做梯度下降
我们有9小时*18项的变量，所以需要9*18个w值，还有一个bias:b
注意下面代码的第三行，这里是整个矩阵前拼接了一列1值，这一列在点乘后就是b
所以这就是为什么dim有18*9+1项
dim = 18 * 9 + 1
w = np.zeros([dim, 1])
x = np.concatenate((np.ones([12 * 471, 1]), x), axis = 1).astype(float)
learning_rate = 100
iter_time = 1000
adagrad = np.zeros([dim, 1])
eps = 0.0000000001
然后定义loss函数，这里知乎作者用的是均方根误差来做loss函数
贴一段百度百科的解释：
均方根误差是预测值与真实值偏差的平方与观测次数n比值的平方根，
在实际测量中，观测次数n总是有限的，真值只能用最可信赖（最佳）值来代替。
标准误差对一组测量中的特大或特小误差反映非常敏感，
所以，标准误差能够很好地反映出测量的精密度。
这正是标准误差在工程测量中广泛被采用的原因。
因此，标准差是用来衡量一组数自身的离散程度，
而均方根误差是用来衡量观测值同真值之间的偏差，
它们的研究对象和研究目的不同，但是计算过程类似。
因为常数项的存在，所以 dimension (dim) 需要多加一栏；
eps 项是避免 adagrad 的分母为 0 而加的极小数值。
rmse公式：
然后adagrid的公式：
这里对gt的计算我没有看懂，根据公式应该是loss对w的偏导，以后搞懂了再来补说明吧
adagrid这一行就是上面的累计求和
然后对w做更新，这样就可以实现梯度下降了
最后保存计算出的w矩阵
for t in range(iter_time):
loss = np.sqrt(np.sum(np.power(np.dot(x, w) - y, 2))/471/12) #rmse
if(t%100==0):
print(str(t) + ":" + str(loss))
gradient = 2 * np.dot(x.transpose(), np.dot(x, w) - y) #dim*1
adagrad += gradient ** 2
w = w - learning_rate * gradient / np.sqrt(adagrad + eps)
np.save('weight.npy', w)
这样我们的w和b就计算完了，接下来拿测试集做测试
还是和train.csv一样先做预处理和Normalize，然后用w和测试集做点乘就可以得到预测值了
记得在拼接1列1作为b
testdata = pd.read_csv('./data/test.csv', header = None, encoding = 'big5')
test_data = testdata.iloc[:, 2:].copy()
test_data.replace('NR', 0, inplace=True)
test_data = test_data.to_numpy()
test_x = np.empty([240, 18*9], dtype = float)
for i in range(240):
test_x[i, :] = test_data[18 * i: 18* (i + 1), :].reshape(1, -1)
for i in range(len(test_x)):
for j in range(len(test_x[0])):
if std_x[j] != 0:
test_x[i][j] = (test_x[i][j] - mean_x[j]) / std_x[j]
test_x = np.concatenate((np.ones([240, 1]), test_x), axis = 1).astype(float)
w = np.load('weight.npy')
ans_y = np.dot(test_x, w)
ans_y[:10]
ans_y就是预测值了，然后按照作业要求写成csv文件。
这里转回pandas处理，因为pandas的to_csv方法很好用
res = pd.DataFrame(ans_y)
res['id'] = ['id_' + str(x) for x in res.index]
res.rename(columns={0:'value'}, inplace=True)
res=res[['id', 'value']] # 调换顺序
res.to_csv('submit.csv', index=False)
最后，可以调整 learning rate、iter_time (iteration 次数)、
取用 features 的多少(取几个小时，取哪些特征)
甚至是不同的 model 来超越 baseline。
