参考自https://blog.csdn.net/qq_46126258/article/details/112468420?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-0&spm=1001.2101.3001.4242

第一次写机器学习的作业，基本上都是借鉴，但是也发现该博主在gradient descent 的 gradient和loss计算时的错误。

实验数据及jupyternotebook 代码

链接：https://pan.baidu.com/s/1koHjoEkpEy7SsLOFRp6Oxg
提取码：ux8i
复制这段内容后打开百度网盘手机App，操作更方便哦

实验环境

• python 3.8
• jupyter notebook + chrome

实验目的

• 利用每天前 9个小时的18个特征，来预测第十个小时的pm2.5
• 给定一个训练集数据和测试集数据，学会对数据的预处理
• 熟悉regression的线性模型，并且用Gradient Descent 和Adamgrad 来优化模型

具体实验步骤

1. 加载数据

# 导入所需要的的包
import sys
import pandas as pd 
import numpy as np
data = pd.read_csv('data/train.csv', encoding = 'big5') 

2. 观察数据
   每天有18项数据，并且RAINFALL存在不能处理的字符串，需要先将字符串更改为可操作的数字，便于操作vector
   

data = data.iloc[:, 3:]  # 提取所有行，并且列从第三列开始
data[data == 'NR'] = 0  # data 和data[] 是不一样的，一个为int，一个为列表
raw_data = data.to_numpy() # 将data

作为一条python新狗，下面连接介绍了pandas 的iloc loc 的具体用法

--------------

3. 特征提取

# 查看 .csv表格，我们观察到，去掉第一行表头，表中有 每月中前二十天的数据，每天24小时记录 18 行feature，
# 所以每天用 18 * 24 矩阵表示
# 一个月有 24 * 20 = 480 个小时，所以12个月由 18 行 * 12 * 480 列

month_data = {} # 按照month做索引，方便以后的数据查找，是一个三维矩阵
for month in range(12):##求12个每月
    sample = np.empty([18, 480])##初始化每月：18 * 480
    for day in range(20):##每个月的20天
        sample[:, day * 24 : (day + 1) * 24] = raw_data[18 * (20 * month + day) : 18 * (20 * month + day + 1), :]
    month_data[month] = sample# 将求好的单月赋值给 month_data
    
    
# 根据题目含义，分为训练集数据和测试集数据，前九个小时做训练集，第十个小时做验证集，所以将十个小时划分为一组。
# 所以可以0-9 为一组，其中0-8为训练集，9为测试集，同时1-10 又为一组
# 因此，可以会得到471组数据，最后一组长9小时，没有验证集 
# 一年有 12 * 471 = 5652 组数据，472 * 12 个测试集，每个测试集有 18 * 9 个feature，对应471 * 12个验证集：PM2.5

x = np.empty([12 * 471, 18 * 9], dtype = float)  # 训练集 12 * 471 行，每列为 18 * 9的数据 
y = np.empty([12 * 471, 1], dtype = float)       #  12 * 471 行， 1列  PM 2.5

#由于已经将一个月的数据合并，那么天与小时都已经是次要的信息。
#因此，此处对于1个月，只需要一个index，从0到471取到就可以了            
for month in range(12):
    for index in range(471):
        x[index + 471 * month, :] = month_data[month][:, index:index + 9].reshape(1, -1)  # 数据reshape为 1 行，列的参数-1会默认得到列
        y[index + 471 * month, :] = month_data[month][9, index + 9]   #  单个数据，一列 ，从第9个位验证数据集，9,10,11,12,13,14 依次为测试集
print(len(x))

4. Normalize 标准化处理数据，就是feature scaling，使得gradient的速度更快

# 正则化处理数据的方法，让数据的值减去均值，再除以标准差
# mean（） 函数求均值，np.std（）求标准差
# 设置的 axis 参数，axis不设置值，就是求所有数的平均值，
# axis = 0， 压缩行，求各列的平均值，返回一个1 * m 的矩阵
# axis = 1，压缩列，求各行的平均值 ，返回一个 n * 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]  # 每个元素减去其所在列的均值，再除以标准差
## 这里就是特征缩放

5. training

# 模型选取：线性模型
# 输入：一次18 * 9个数据
# 参数：18 * 9个weight + bias
# loss function：均方根误差
# 优化算法：gradient descent，Adagrad 
# np.concatenate((a,b),axis=1) 进行拼接，拼接结果是[a,b], 相当于在x矩阵的左侧第一行添加同等行数的1
# np.dot(a,b) a与b进行矩阵乘法 
import matplotlib.pyplot as plt
dim = 18 * 9 + 1 # w的行数，w为 18 * 9 + 1 行
w = np.ones([dim, 1]) 
x = np.concatenate( (np.ones([12 * 471, 1]), x), axis=1).astype(float)  # 每运行一次加一列
print(x.shape)
# 对列进行拼接，x矩阵被拼接为 471 * 12行，18*9 + 1列，加上的一列为bias
#x = np.concatenate( (np.ones([12 * 471, 1]), x), axis=1).astype(float)
N = 12 * 471 # y-hat 的行的行数
learning_rate = 0.03 
iter_time = 1000 # 训练的迭代次数
loss = np.zeros([iter_time, 1])  # 每次迭代loss更新，所以与迭代次数相同
for i in range(iter_time):
    y1 = np.dot(x, w) # x 为 471*12行，18*9 +1 列 ，w为 18*9+1行 1列
    #loss[i] = np.sqrt(np.sum(np.power(y1-y,2)/N))  
    #gradient = (2/N) * np.dot(x.T,y1-y)  # dot 为两个矩阵相乘 x.t (163,5652) * (5652, 1) = (163,1)
    # 上面的loss function 和 gradient 的计算是不匹配的，下面的匹配 
    loss[i] = np.sum(np.power(y1-y, 2)/N )
    gradient = (2/N) * np.dot(x.T,y1-y)
    w = w - learning_rate * gradient 
plt.plot(loss)  # 为了visualize loss，来调整learning_rate
plt.show()

6. AdamGrad 优化

## AdaGrad优化
## 上一个算法中learning_rate是凭借试验和经验设置的，且是一个常数保持不变，所以在计算loss的时候可能卡在一个曲线左右循坏跳，而到不了minLoss
## 那么可不可以让他自己根据每一步求出的w来进行自己的更新，自动调节learning_rate
## AdaGrad算法，一开始是激励收敛，到后面变为惩罚收敛，速度越来越慢
x = np.delete(x, 0, axis = 1)
print(x[0][0])
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
for t in range(iter_time):
    loss = np.sqrt(np.sum(np.power(np.dot(x, w) - y,2)) / 471 / 12)
    if(t%100==0):
        print(str(t) + ":" + str(loss))
    gradient = 2 * np.dot(x.transpose(), np.dot(x,w) - y)
    adagrad += gradient ** 2 #求adagrad的平方power(adagrad,2)  用一次微分的和来替代二次微分（见下图）
    w = w - learning_rate * gradient / np.sqrt(adagrad + eps)
np.save('weight.npy',w)

7. testing

test_data = pd.read_csv('data/test.csv', header = None, encoding = 'big5')
test_data = test_data.iloc[:, 2:] # 取第二列到最后一列的所以行 
test_data[test_data == 'NR'] = 0
test_data = test_data.to_numpy() # 转换为numpy数据
test_x = np.empty([240, 18*9],dtype=float)
for i in range(240):#将每天的18项数据横排摆放
    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:#确保分母不为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) # 增加一列数据计算bias权重

w = np.load('weight.npy')
ans_y = np.dot(test_x, w)
ans_y
import csv
with open('submit.csv', mode='w', newline='') as submit_file:
    csv_writer = csv.writer(submit_file)
    header = ['id', 'value']
    print(header)
    csv_writer.writerow(header)
    for i in range(240):
        row = ['id_' + str(i), ans_y[i][0]]
        csv_writer.writerow(row)
        print(row)