Hung-Yi Lee homework[1]:regression

---

  鸣谢：李宏毅2020机器学习作业1——Linear Regression

  train.csv,test.csv提取码：42fs

  全部代码如下：

import sys
import pandas as pd
import numpy as np
import math

# 导入数据
data = pd.read_csv('train.csv', encoding='big5')

# 分离出第四列开始的数据
data = data.iloc[:, 3:]
data[data == 'NR'] = 0
raw_data = data.to_numpy()
print('raw_data')
print(raw_data)

# 对data进行调整，将4320*24重组为12*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*(20*month+day):18*(20*month+day+1), :]
        month_data[month] = sample

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]
print('x:')
print(x)
print('y:')
print(y)

# 归一化
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]
# 将训练集分成训练-验证集，用来最后检验我们的模型
x_train_set = x[: math.floor(len(x) * 0.8), :]
y_train_set = y[: math.floor(len(y) * 0.8), :]
x_validation = x[math.floor(len(x) * 0.8):, :]
y_validation = y[math.floor(len(y) * 0.8):, :]
print(x_train_set)
print(y_train_set)
print(x_validation)
print(y_validation)
print(len(x_train_set))
print(len(y_train_set))
print(len(x_validation))
print(len(y_validation))

# 因为存在偏差bias，所以dim+1
dim = 18 * 9 + 1
# w维度为163*1
w = np.zeros([dim, 1])
# x_train_set维度为 4521*163
x_train_set = np.concatenate((np.ones([len(x_train_set), 1]), x_train_set), axis=1).astype(float)
# 设置学习率
learning_rate = 10
# 设置迭代数
iter_time = 30000
# RMSprop参数初始化
adagrad = np.zeros([dim, 1])
eps = 0.0000000001
# beta = 0.9
# 迭代
for t in range(iter_time):
    loss = np.sqrt(np.sum(np.power(np.dot(x_train_set, w)-y_train_set, 2))/len(x_train_set))
    if t % 100 == 0:
        print("迭代的次数：%i ， 损失值：%f" % (t, loss))
        # gradient = 2*np.dot(x.transpose(),np.dot(x,w)-y)
    # 计算梯度值
    gradient = (np.dot(x_train_set.transpose(), np.dot(x_train_set, w)-y_train_set))/(loss*len(x_train_set))
    adagrad += (gradient ** 2)
    # 更新参数w
    w = w - learning_rate * gradient / np.sqrt(adagrad + eps)
# 保存参数w
np.save('weight.npy', w)

# 读入测试数据集并对其进行预处理
testdata = pd.read_csv('test.csv', header=None, encoding='big5')
test_data = testdata.iloc[:, 2:]
test_data[test_data == 'NR'] = 0
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)
print(test_x)

# 验证模型并预测
x_validation = np.concatenate((np.ones([len(x_validation), 1]), x_validation), axis=1).astype(float)
for m in range(len(x_validation)):
    Loss = np.sqrt(np.sum(np.power(np.dot(x_validation, w)-y_validation, 2))/len(x_validation))
print("the loss on val data is %f" % Loss)

# 预测
ans_y = np.dot(test_x, w)
print("预测P,2.5值")
print(ans_y)

  开始学习了才发现，自己对于python一无所知，其实之前有断断续续写过一点python的代码，但都没有登堂入室，这次也算是用到啥学一点，各位就不要嫌弃这篇文章啰嗦了。

1. 导入所需要的库

import sys
import pandas as pd
import numpy as np
import math

  在后续的实验中，可以发现我们并没有使用sys库，所以第一句import sys也是可以不用的。

  在pycharm中下载安装numpy后，执行代码的时候遇到了一个runtime error【因为这篇博客是在跑通代码后的总结和思考，所以具体错误并没有截图，在历史记录中找了遇到相同问题博主的截图】

RuntimeError: The current Numpy installation ('D:\\ProgramData\\Python\\Python39\\lib\\site-packages\\numpy\\__init__.py') fails to pass a sanity check due to a bug in the windows runtime. See this issue for more information: https://tinyurl.com/y3dm3h86

  解决方案是：将numpy卸载安装一个老版本(我后来使用的版本：numpy==1.19.3)

numpy库：python的一个扩展程序库，支持大量的维度数组与矩阵运算

pandas库：一个强大的分析结构化数据的工具集

math库：数学运算的库

2. 导入train数据

# 导入数据
data = pd.read_csv('train.csv', encoding='big5')

  通过对数据文件的分析，train.csv的资料为12个月，每个月取20天，每天24小时的资料（每个小时的资料有18个features）

  所以我们需要把无用信息都剔除出去

  对数据进行处理，取第4列开始的数据【从这里我们可以看出，从pandas的read_csv()函数中读取出的数据是从0开始计数的】

# 分离出第四列开始的数据
data = data.iloc[:, 3:]
data[data == 'NR'] = 0
raw_data = data.to_numpy()
print('raw_data')
print(raw_data)

pandas.DataFrame.iloc

  .iloc[] is primarily integer position based (from 0 to length-1 of the axis), but may also be used with a boolean array.

  Allowed inputs are:

• An integer 5
• A list or Array of integers [4,3,0]
• A slice object with ints 1:7
• A boolean array
• A callable function with one argument(the calling Series or DataFrame)and that returns valid output for indexing (one of the above). This is useful in method chains, when you don’t have a reference to the calling object, but would like to base your selection on some value

  .iloc如果请求的索引器超出范围，.iloc将引发indexError,但切片索引器除外，该索引器允许超出范围的索引。

pandas.DataFrame.to_numpy

  DataFrame.to_numpy(dtype=None,copy=False,na_value=)

  作用：convert the DataFrame to a NumPy array.

  na_value:用于missing values。

Numpy数组

  Numpy中定义的最重要的对象是称为ndarray的N维数组类型，它是描述相同类型的元素集合。ndarray中的每个元素都是数据类型对象(dtype)的对象。ndarray中的每个元素在内存中使用相同大小的块。

numpy.array(object, dtype=None, copy=True, order='K', subok=False, ndmin=0)

  ndmin:数组的最小维数。

  可以看到结果已经变成了二维数组。

ndarray的内部机理

  ndarray的组成：

  ndarray与数组不同，它不仅仅包含数据信息，还包括其他描述信息。ndarray内部由以下内容组成：

• 数据指针：一个指向实际数据的指针。
• 数据类型(dtype)：描述了每个元素所占字节数。
• 维度(shape)：一个表示数组形状的元组。
• 跨度(stride)：一个表示从当前维度前进到下一维度的当前位置所需要"跨过"的字节数。

  numpy中，数据存储在一个均匀连续的内存块中，可以这么理解，numpy将多维数组在内部以一维数组的方式存储，我们只要知道了每个元素所占的字节数(dtype)以及每个维度中元素的个数(shape)，就可以快速定位到任意维度的任意一个元素。

3. data数据的预处理

  我们需要对数据进行重组，以便进行我们的训练。

# 对data进行调整，将4320*24重组为12*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*(20*month+day):18*(20*month+day+1), :]
        month_data[month] = sample

  按照作业要求，每9个小时的数据来预测第10个小时的PM2.5，每天24个小时，每9个小时构成一个data，第10个小时为label，每天有$24-10+1=15$个data和label，每个月有300个data，所以一年有$12\times 20\times 15=3600$个data

  因为一个月有20天都是连续的，所以可以将20天的280个小时都看成是连续的，所以一个月有$480-10+1=471$个data，一年有$471\times 12=5652$个data,同样有5652个label。采用这种方法我们可以构建比较多的data。

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]

  其中有一句if day == 19 and hour>14这句话代表的是某一个月的第20天的最后一组数据了（因为月和月之间不是连续的）。

reshape(1,-1)

In [2]: arr=np.arange(16).reshape(2,8)
out[2]:
 
In [3]: arr
out[3]:
array([[ 0,  1,  2,  3,  4,  5,  6,  7],
       [ 8,  9, 10, 11, 12, 13, 14, 15]])
 
In [4]: arr.reshape(4,-1) #将arr变成4行的格式，列数自动计算的(c=4, d=16/4=4)
out[4]:
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])
 
In [5]: arr.reshape(8,-1) #将arr变成8行的格式，列数自动计算的(c=8, d=16/8=2)
out[5]:
array([[ 0,  1],
       [ 2,  3],
       [ 4,  5],
       [ 6,  7],
       [ 8,  9],
       [10, 11],
       [12, 13],
       [14, 15]])
 
In [6]: arr.reshape(10,-1) #将arr变成10行的格式，列数自动计算的(c=10, d=16/10=1.6 != Int)
out[6]:
ValueError: cannot reshape array of sie 16 into shape (10,newaxis)

4. normalization

# 归一化
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]

  归一化是为了进行feature scaling，feature scaling使得参数的更新更为快速且简单。

5. 从训练集中取出一部分设立验证集

# 将训练集分成训练-验证集，用来最后检验我们的模型
x_train_set = x[: math.floor(len(x) * 0.8), :]
y_train_set = y[: math.floor(len(y) * 0.8), :]
x_validation = x[math.floor(len(x) * 0.8):, :]
y_validation = y[math.floor(len(y) * 0.8):, :]

  从训练集中取出了80%作为训练集。

6. 开始训练

# 因为存在偏差bias，所以dim+1
dim = 18 * 9 + 1
# w维度为163*1
w = np.zeros([dim, 1])
# x_train_set维度为 4521*163
x_train_set = np.concatenate((np.ones([len(x_train_set), 1]), x_train_set), axis=1).astype(float)
# 设置学习率
learning_rate = 10
# 设置迭代数
iter_time = 30000
# RMSprop参数初始化
adagrad = np.zeros([dim, 1])
eps = 0.0000000001
# beta = 0.9
# 迭代
for t in range(iter_time):
    loss = np.sqrt(np.sum(np.power(np.dot(x_train_set, w)-y_train_set, 2))/len(x_train_set))
    if t % 100 == 0:
        print("迭代的次数：%i ， 损失值：%f" % (t, loss))
        # gradient = 2*np.dot(x.transpose(),np.dot(x,w)-y)
    # 计算梯度值
    gradient = (np.dot(x_train_set.transpose(), np.dot(x_train_set, w)-y_train_set))/(loss*len(x_train_set))
    adagrad += (gradient ** 2)
    # 更新参数w
    w = w - learning_rate * gradient / np.sqrt(adagrad + eps)
# 保存参数w
np.save('weight.npy', w)

• 设置超参数：学习率、迭代次数登
• 计算损失L
• 计算梯度gradient
• 梯度下降

  损失函数采用均方根误差公式：
$$L=\sqrt{\frac{\sum (xw-y{)}^2}{num}}$$
  对参数$w$计算梯度值：
$$dw=\frac{x(xw-y)}{L\times num}$$

adagrad(Adaptive gradient alorithm)

  这个算法可以对低频的参数做交大的更新，对高频的参数做较小的更新，也因此，他对系数的数据表现很好，很好地提高了SGD的鲁棒性。

  权值更新的规则：
$${\theta }_{t+1,i}={\theta }_{t,i}-\frac{\eta }{\sqrt{G_{t,ii}+ϵ}}\cdot g_{t,i}$$
  其中，$g$为：t时刻${\theta }_i$的梯度。$G$是一个对角矩阵，$(i,i)$元素就是$t$时刻参数${\theta }_i$前面所有时刻的梯度平方和。

  Adagrad 的优点是减少了学习率的手动调节。

  Adagrad的缺点是分母会不断积累，这样学习率就会收缩并最终会变得非常小。

7. 对测试集进行预处理

# 读入测试数据集并对其进行预处理
testdata = pd.read_csv('test.csv', header=None, encoding='big5')
test_data = testdata.iloc[:, 2:]
test_data[test_data == 'NR'] = 0
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)

8. 验证模型并预测

# 验证模型并预测
x_validation = np.concatenate((np.ones([len(x_validation), 1]), x_validation), axis=1).astype(float)
for m in range(len(x_validation)):
    Loss = np.sqrt(np.sum(np.power(np.dot(x_validation, w)-y_validation, 2))/len(x_validation))
print("the loss on val data is %f" % Loss)

# 预测
ans_y = np.dot(test_x, w)
print("预测P,2.5值")
print(ans_y)

  【要点】在看到代码的时候，不经拍案叫绝，采用了验证集的loss来代表测试集的loss。

  这是因为验证集和训练集本身也毫不相关，在和训练集的关系上，验证集和测试集处于同等的地位，所以用验证集的loss来代替测试集的loss是可行的。