线性回归梯度下降求解—代码实现

李毅宏机器学习——课后作业Linear Regression

使用某个气象站的观测记录，分为train set训练集和test set测试集。训练集是一年中每个月的前20天每天24小时的所有资料，测试集是从剩下的资料中取样出来。
Data中含有18项观测数据包括PM2.5、PM10等，模型的目标是通过前9个小时的观测数据，预测第10个小时的PM2.5指标，显然Linear Regression线性回归模型最适合。
Regression回归的理论基础，Regression回归这篇笔记中已经详细讨论过。

作业的数据在链接：https://pan.baidu.com/s/1_26vBKu3RHSD9YfJZBP67A
提取码：fv58

一、数据预处理

数据预处理部分代码：

import sys
import pandas as pd 
import numpy as np 
data=pd.read_csv('./train.csv',encoding='big5')
data=data.iloc[:,3:]
data[data=='NR']=0#将'NR'都替换成0
raw_data=np.array(data.values.tolist())#将数据转成二维数组

将数据从纵向时序排列调整成横向18个特征*一年12个月20天每天24小时的数据组。

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

在横向时间轴上以每10小时为一个训练集滑动提取训练集数据

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) #vector dim:18*9 (9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9)
            y[month * 471 + day * 24 + hour, 0] = month_data[month][9, day * 24 + hour + 9] #value

到这里数据集就转化为：

二、特征缩放

数据处理好后不能直接跑模型，可以先进行特征归一化平衡特征值，能够使得梯度更加容易得下降到全局最优。

其中mi是特征的平均值，σi是特征最大值减去最小值。

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]

将训练集按8：2比例分为组内训练集与测试集，测试模型训练效果。

import math
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): , :]

三、梯度下降求解

梯度下降法是线性回归模型参数求解的方法，在Gradient descent 梯度下降中讨论了计算方法与优化方法。

dim = 18 * 9 + 1 #之所以+1是留出误差项b的位置
w = np.ones([dim, 1]) #w就是待学习参数
x = np.concatenate((np.ones([12 * 471, 1]), x), axis = 1).astype(float)
learning_rate = 100 #初始化学习率设置为100
iter_time = 1000 #训练次数为1000次
adagrad = np.zeros([dim, 1]) #adagrad自适应学习率方法
eps = 0.0000000001 #adagrad中防止分母为0的极小值
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):#每隔100次输出一次误差值
        print(str(t) + ":" + str(loss))
    gradient = 2 * np.dot(x.transpose(), np.dot(x, w) - y) #计算的梯度，在下文中详解
    adagrad += gradient ** 2 #adagrad中，学习率要除以历史所有的梯度平方和
    w = w - learning_rate * gradient / np.sqrt(adagrad + eps) #更新w的值
np.save('weight.npy', w) #储存w

这里要详细解释为什么代码中梯度gradient的计算公式是：

gradient = 2 * np.dot(x.transpose(), np.dot(x, w) - y)

即：

首先，线性回归的表达式为：

令X0=1，则表达式可以化为：

损失函数表达式：

对w0求微分：

这样求完所有参数w的微分，就是：

至此可以看出，梯度gradient就是2×X的转置×损失误差。

至于adagrad的表达式：

模型训练过程输出可以看到损失误差逐步减小。

四、测试模型效果

w = np.load('weight.npy') #提取训练好的参数w
x_validation = np.concatenate((np.ones([1131, 1]), x_validation), axis=1).astype(float)
ans_y = np.dot(x_validation, w)
loss = np.sqrt(np.sum(np.power(ans_y - y_validation, 2)) / 1131)
print(loss) #输出误差

测试数据的误差为6.39619，是个还不错的结果。

五、提交作业

用同样的预处理方法处理test set测试集数据，输出预测结果保存到csv文件。

testdata = pd.read_csv('./test.csv', header=None, encoding='big5')
test_data = testdata.iloc[:, 2:]  
test_data[test_data == 'NR'] = 0  
test_data = np.array(test_data.values.tolist())  
test_x = np.empty([240, 18 * 9], dtype=float)  
for i in range(240):  # 共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)  # test data的预测值ans_y=test_x与W的积
#输出到csv文件
import csv
with open('F:/BaiduNetdiskDownload/hw1/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)

完整做完这套作业，对于不理解的细节进行推敲，感觉自己对于逻辑回归与梯度下降有了更多实质的理解。