机器学习与数据挖掘--梯度下降法训练多元线性回归模型

机器学习与数据挖掘实验二

（使用梯度下降法训练多元线性回归模型）

实验目的：

熟练掌握线性回归的基本原理以及梯度下降法和最小二乘法

实验环境：

Anaconda/Jupyter notebook/Pycharm

实验内容：

1. 编码实现基于梯度下降的多变量线性回归算法，包括梯度的计算与验证；
2. 画数据散点图，以及得到的直线；
3. 画梯度下降过程中损失的变化图；
4. 基于训练得到的参数，输入新的样本数据，输出预测值；

实验步骤：

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
df_ads = np.loadtxt('./data2.txt',delimiter=',')

x = np.array(df_ads)
y = np.array(df_ads[:,2])
print("张量x的阶：",x.ndim)
print("张量x的形状",x.shape)
print(x)

y = y.reshape(-1,1)
print(y.shape)

from sklearn import datasets
from mpl_toolkits.mplot3d import Axes3D
x1 = x[:,0]#自变量1
x2 = x[:,1]#自变量2
x3 = x[:,-1]#因变量
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x1,x2,x3)
plt.show()

# 将数据继进行80%（训练集）和20%（验证集）的分割
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y,test_size=0.2, random_state=0)

def scaler(train,test): #  定义归一化函数。进行数据压缩
        #数据的压缩
        min = train.min(axis=0) # 训练集最小值
        max = train.max(axis=0) # 训练集最大值
        gap = max - min # 最大值和最小值的差
        train -= min # 所有数据减最小值
        train /= gap # 所有数据除以大小值差
        test -= min #把训练集最小值应用于测试集
        test /= gap #把训练集大小值差应用于测试集
        return train, test # 返回压缩后的数据

def min_max_gap(train): # 计算训练集最大，最小值以及他们的差，用于后面反归一化过程
    min = train.min(axis=0) # 训练集最小值
    max = train.max(axis=0) # 训练集最大值
    gap = max - min # 最大值和最小值的差
    return min, max, gap   
y_min, y_max, y_gap = min_max_gap(y_train)

x_train_original = x_train.copy() # 保留一份训练集数据副本，用于对要预测数据归一化

x_train,x_test = scaler(x_train,x_test) # 对特征归一化
y_train,y_test = scaler(y_train,y_test) # 对标签也归一化

x0_train = np.ones((len(x_train),1)) # 构造X_train长度的全1数组配合对Bias的点积
x_train = np.append(x0_train, x_train, axis=1) #把X增加一系列的1
x0_test = np.ones((len(x_test),1)) # 构造X_test长度的全1数组配合对Bias的点积
x_test = np.append(x0_test, x_test, axis=1) #把X增加一系列的1
print ("张量X的形状:", x_train.shape)
print (x_train)

def loss_function(x, y, W): # 手工定义一个MSE均方误差函数,W此时是一个向量
    y_hat = x.dot(W.T) # 点积运算 h(x)=w_0*x_0 + w_1*x_1 + w_2*x_2 + w_3*x_3    
    print("x.shape",x.shape)
    print("W.shape",W.shape)
    loss = y_hat.reshape((len(y_hat),1))-y # 中间过程,求出当前W和真值的差异
    cost = np.sum(loss**2)/(2*len(x)) # 这是平方求和过程, 均方误差函数的代码实现
    return cost # 返回当前模型的均方误差值

def gradient_descent(x, y, W, lr, iterations): # 定义梯度下降函数
    l_history = np.zeros(iterations) # 初始化记录梯度下降过程中损失的数组
    W_history = np.zeros((iterations,len(W))) # 初始化权重数组 
    for iter in range(iterations): # 进行梯度下降的迭代，就是下多少级台阶
        y_hat = x.dot(W) # 这个是向量化运行实现的假设函数   
        loss = y_hat.reshape((len(y_hat),1))-y # 中间过程, y_hat和y真值的差
        derivative_W = x.T.dot(loss)/(2*len(x)) #求出多项式的梯度向量
        derivative_W = derivative_W.reshape(len(W)) 
        W = W - alpha*derivative_W # 结合下降速率更新权重
        l_history[iter] = loss_function(x, y, W) # 损失的历史记录 
        W_history[iter] = W # 梯度下降过程中权重的历史记录
    return l_history, W_history # 返回梯度下降过程数据

#首先确定参数的初始值
iterations = 300; # 迭代300次
alpha = 0.15; #学习速率设为0.15
weight = np.array([0.5,1,1,1]) # 权重向量，w[0] = bias
#计算一下初始值的损失
print ('当前损失：',loss_function(x_train, y_train, weight)

# 定义线性回归模型
def linear_regression(x, y, weight, alpha, iterations): 
    loss_history, weight_history = gradient_descent(x, y, 
                                                    weight, 
                                                    alpha, iterations)
    print("训练最终损失:", loss_history[-1]) # 打印最终损失
    y_pred = x.dot(weight_history[-1]) # 进行预测
    traning_acc = 100 - np.mean(np.abs(y_pred - y))*100 # 计算准确率
    print("线性回归训练准确率: {:.2f}%".format(traning_acc))  # 打印准确率
    return loss_history, weight_history # 返回训练历史记录

# 调用刚才定义的线性回归模型
loss_history, weight_history = linear_regression(x_train, y_train,weight, alpha, iterations) #训练机器

print("权重历史记录：", weight_history)
print("损失历史记录：", loss_history)

plt.plot(loss_history,'g--',label='Loss Curve')
plt.xlabel('Iterations') # x轴Label
plt.ylabel('Loss') # y轴Label
plt.legend() # 显示图例
plt.show() # 显示损失曲线
print(alpha)

x_plan = [250,50,50] # 要预测的X特征数据
x_train,x_plan = scaler(x_train_original,x_plan) # 对预测数据也要归一化缩放
x_plan = np.append([1], x_plan ) # 加一个哑特征X0 = 1
y_plan = np.dot(weight_history[-1],x_plan) # [-1] 即模型收敛时的权重
# 对预测结果要做反向缩放，才能得到与原始广告费用对应的预测值
y_value = y_plan*y_gap + y_min # y_gap是当前y_train中最大值和最小值的差，y_min是最小值
print ("预计商品销售额： ",y_value, "千元") 

x_p=x.copy()
x_t_o=x_train_original.copy()
x_t,x_p=scaler(x_t_o,x_p)
c=np.ones(x_p.shape[0]).transpose()
x_p=np.insert(x_p,0,values=c,axis=1)
y_p=x_p.dot(weight_history[-1])
y_v=y_p*y_gap+y_min

fig=plt.figure()
ax=Axes3D(fig)

x_p0=np.array(x_p[:,0])
x_p1,x_p2=np.meshgrid(x_p[:,1],x_p[:,2])
y_p1=weight_history[-1][0]+x_p1.dot(weight_history[-1][1])+x_p2.dot(weight_history[-1][2])
y_v1=y_p1*y_gap+y_min

x_p1,x_p2=np.meshgrid(x[:,0],x[:,1])
ax.plot_surface(x_p1,x_p2,y_v1,cmap='hot')