线性回归，梯度下降

最新推荐文章于 2022-05-01 18:22:33 发布

氰镁JUNK

最新推荐文章于 2022-05-01 18:22:33 发布

阅读量258

点赞数 1

文章标签：推荐系统

本文链接：https://blog.csdn.net/L2350519620/article/details/115792141

版权

线性回归

1.线性回归简介

线性回归是利用数理统计中的回归分析，来确定两种或两种以上变量间相互依赖的变量关系的一种统计分析方法，运用十分广泛。
回归分析中，如果只包括一个自变量和一个因变量，且两者的关系可以用一条直线来表示，这种回归分析成为一元线性回归分析；如果回归分析中包括两个或两个以上的变量，且因变量和自变量之间是先线性关系，则称为多元线性回归分析。

2.线性回归分析的一般形式

在这里插入图片描述

3.回归分析可以解决的问题

1.建立变量之间的数学公式，通常称为经验公式。
2.利用概率统计基础知识进行分析，从而判断所建立的经验公式的有效性。
3.进行因素分析，确定影响某一变量的若干因素中，何者为主要，何者为次要，以及它们之间的关系

一元线性回归

1，函数的推导过程

设线性回归函数：在这里插入图片描述
构造损失函数：
思路：通过梯度下降不断更新θ0和θ1，当损失函数的值特别小的时候，就得到了我们的函数模型。
过程：

step1.求导

step2.更新θ0和θ1 在这里插入图片描述
step3.带入损失函数，求损失函数的值，若得到的值小于 $\xi$ (一般为0.01或者0.001这样的小数），退出；否则，返回step1。

2，代码实现

批量梯度下降法

import matplotlib.pyplot as plt
import matplotlib
from math import pow
from random import uniform

x = [1,3,5,7,9,11,13]
y = [100,111,130,144,149,166,179]

#线性回归函数为 y=theta0+theta1*x
#参数定义
theta0 = uniform(0,2)#对 theata0 随机赋值
theta1 = uniform(0,2)#对 theata1 随机赋值
alpha = 0.1#学习率
m = len(x)
count = 0
loss = []

for num in range(10000):
    count += 1
    diss = 0   #误差
    deriv0 = 0 #对 theata0 导数
    deriv1 = 0 #对 theata1 导数
    #求导
    for i in range(m):
        deriv0 += (theta0+theta1*x[i]-y[i])/m
        deriv1 += ((theta0+theta1*x[i]-y[i])/m)*x[i]
    #更新 theta0 和 theta1
    for i in range(m):
        theta0 = theta0 - alpha*((theta0+theta1*x[i]-y[i])/m) 
        theta1 = theta1 - alpha*((theta0+theta1*x[i]-y[i])/m)*x[i]
    #求损失函数 J (θ)
    for i in range(m):
        diss = diss + (1/(2*m))*pow((theta0+theta1*x[i]-y[i]),2)
    loss.append(diss)

    #如果误差已经很小，则退出循环
    if diss <= 0.001:
        break
    
print("本次迭代次数为：{}次，最终得到theta0={}，theta1={}".format(count,theta0,theta1))
print("本次迭代得到的回归函数是：y={}+{}*x".format(theta0,theta1))
#画原始数据图和得到的线性回归函数图
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
plt.plot(x,y,'bo',label='数据')
plt.plot(x,[theta0+theta1*x for x in x],label='线性回归函数')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
#画损失函数（误差）变化图
plt.scatter(range(count),loss)
plt.show()

输出结果：
本次迭代次数为：10000次，最终得到theta0=96.86808644268262，theta1=5.762687172041142
本次迭代得到的回归函数是：y=96.86808644268262+5.762687172041142*x
在这里插入图片描述

2，随机梯度先将

import matplotlib.pyplot as plt
import matplotlib
from math import pow
import random

x = [1,3,5,7,9,11,13]
y = [100,111,130,144,149,166,179]

#线性回归函数为 y=theta0+theta1*x
#参数定义
theta0 = random.uniform(0,2)#对 theata0 随机赋值
theta1 = random.uniform(0,2)#对 theata1 随机赋值
alpha = 0.1#学习率
m = len(x)
count = 0
loss = []

for num in range(10000):
    count += 1
    diss = 0   #误差
    deriv0 = 0 #对 theata0 导数
    deriv1 = 0 #对 theata1 导数
    #求导
    for i in range(m):
        deriv0 += (theta0+theta1*x[i]-y[i])/m
        deriv1 += ((theta0+theta1*x[i]-y[i])/m)*x[i]
    #更新 theta0 和 theta1
    for i in range(m):
        theta0 = theta0 - alpha*((theta0+theta1*x[i]-y[i])/m) 
        theta1 = theta1 - alpha*((theta0+theta1*x[i]-y[i])/m)*x[i]
    #求损失函数 J (θ)：
    rand_i = random.randrange(0,m)
    diss = diss + (1/(2*m))*pow((theta0+theta1*x[rand_i]-y[rand_i]),2)
    loss.append(diss)

    #如果误差已经很小，则退出循环
    if diss <= 0.001:
        break
    
print("本次迭代次数为：{}次，最终得到theta0={}，theta1={}".format(count,theta0,theta1))
print("本次迭代得到的回归函数是：y={}+{}*x".format(theta0,theta1))
#画原始数据图和得到的线性回归函数图
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
plt.plot(x,y,'bo',label='数据')
plt.plot(x,[theta0+theta1*x for x in x],label='线性回归函数')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
#画损失函数（误差）变化图
plt.scatter(range(count),loss)
plt.show()

输出结果为：
本次迭代次数为：10000次，最终得到theta0=96.86808644268262，theta1=5.762687172041142
本次迭代得到的回归函数是：y=96.86808644268262+5.762687172041142*x
在这里插入图片描述

注意：因为是随机的，所以每次结果都不一样。
3，小批量梯度下降

import matplotlib.pyplot as plt
import matplotlib
from math import pow
import random

x = [1,3,5,7,9,11,13]
y = [100,111,130,144,149,166,179]

#目标函数为 y=theta0+theta1*x
#参数定义
theta0 = random.uniform(0,2)#对 theata0 随机赋值
theta1 = random.uniform(0,2)#对 theata1 随机赋值
alpha = 0.1#学习率
m = len(x)
count = 0
loss = []

for num in range(10000):
    count += 1
    diss = 0   #误差
    deriv0 = 0 #对 theta0 导数
    deriv1 = 0 #对 theta1 导数
    #求导
    for i in range(m):
        deriv0 += (theta0+theta1*x[i]-y[i])/m
        deriv1 += ((theta0+theta1*x[i]-y[i])/m)*x[i]
    #更新 theta0 和 theta1
    for i in range(m):
        theta0 = theta0 - alpha*((theta0+theta1*x[i]-y[i])/m) 
        theta1 = theta1 - alpha*((theta0+theta1*x[i]-y[i])/m)*x[i]
    #求损失函数 J (θ)
    rand_ls = list(range(3))
    for i in range(3):
        rand_i = random.randrange(0,m)
        rand_ls[i] = rand_i
    for i in rand_ls:
        diss = diss + (1/(2*m))*pow((theta0+theta1*x[i]-y[i]),2)
    loss.append(diss)

    #如果误差已经很小，则退出循环
    if diss <= 0.001:
        break
    
print("本次迭代次数为：{}次，最终得到theta0={}，theta1={}".format(count,theta0,theta1))
print("本次迭代得到的回归函数是：y={}+{}*x".format(theta0,theta1))
#画原始数据图和得到的线性回归函数图
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
plt.plot(x,y,'bo',label='数据')
plt.plot(x,[theta0+theta1*x for x in x],label='线性回归函数')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
#画损失函数（误差）变化图
plt.scatter(range(count),loss)
plt.show()

输出结果为：
本次迭代次数为：10000次，最终得到theta0=96.86808644268262，theta1=5.762687172041142
本次迭代得到的回归函数是：y=96.86808644268262+5.762687172041142*x
在这里插入图片描述

在这里插入图片描述

3，举例

波士顿房价预测
在这里插入图片描述

#房屋价格与面积
#序号：1     2    3    4     5     6     7
#面积：150  200  250  300   350   400   600  
#价格：6450 7450 8450 9450 11450 15450 18450 

import matplotlib.pyplot as plt
import matplotlib
from math import pow
from random import uniform
import random

x0 = [150,200,250,300,350,400,600]
y0 = [6450,7450,8450,9450,11450,15450,18450]
#为了方便计算，将所有数据缩小 100 倍
x = [1.50,2.00,2.50,3.00,3.50,4.00,6.00]
y = [64.50,74.50,84.50,94.50,114.50,154.50,184.50]


#线性回归函数为 y=theta0+theta1*x
#参数定义
theta0 = 0.1#对 theata0 赋值
theta1 = 0.1#对 theata1 赋值
alpha = 0.1#学习率
m = len(x)
count0 = 0
theta0_list = []
theta1_list = []

#使用批量梯度下降法
for num in range(10000):
    count0 += 1
    diss = 0   #误差
    deriv0 = 0 #对 theata0 导数
    deriv1 = 0 #对 theata1 导数
    #求导
    for i in range(m):
        deriv0 += (theta0+theta1*x[i]-y[i])/m
        deriv1 += ((theta0+theta1*x[i]-y[i])/m)*x[i]
    #更新 theta0 和 theta1
    for i in range(m):
        theta0 = theta0 - alpha*((theta0+theta1*x[i]-y[i])/m) 
        theta1 = theta1 - alpha*((theta0+theta1*x[i]-y[i])/m)*x[i]
    #求损失函数 J (θ)
    for i in range(m):
        diss = diss + (1/(2*m))*pow((theta0+theta1*x[i]-y[i]),2)
        
    theta0_list.append(theta0*100)
    theta1_list.append(theta1)
    #如果误差已经很小，则退出循环
    if diss <= 0.001:
        break
    
theta0 = theta0*100#前面所有数据缩小了 100 倍，所以求出的 theta0 需要放大 100 倍，theta1 不用变

#使用随机梯度下降法
theta2 = 0.1#对 theata2 赋值
theta3 = 0.1#对 theata3 赋值
count1 = 0
theta2_list = []
theta3_list = []

for num in range(10000):
    count1 += 1
    diss = 0   #误差
    deriv2 = 0 #对 theata2 导数
    deriv3 = 0 #对 theata3 导数
    #求导
    for i in range(m):
        deriv2 += (theta2+theta3*x[i]-y[i])/m
        deriv3 += ((theta2+theta3*x[i]-y[i])/m)*x[i]
    #更新 theta0 和 theta1
    for i in range(m):
        theta2 = theta2 - alpha*((theta2+theta3*x[i]-y[i])/m) 
        theta3 = theta3 - alpha*((theta2+theta3*x[i]-y[i])/m)*x[i]
    #求损失函数 J (θ)
    rand_i = random.randrange(0,m)
    diss = diss + (1/(2*m))*pow((theta2+theta3*x[rand_i]-y[rand_i]),2)

    theta2_list.append(theta2*100)
    theta3_list.append(theta3)
    #如果误差已经很小，则退出循环
    if diss <= 0.001:
        break
theta2 = theta2*100
    
print("批量梯度下降最终得到theta0={}，theta1={}".format(theta0,theta1))
print("           得到的回归函数是：y={}+{}*x".format(theta0,theta1))
print("随机梯度下降最终得到theta0={}，theta1={}".format(theta2,theta3))
print("           得到的回归函数是：y={}+{}*x".format(theta2,theta3))
#画原始数据图和函数图
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
plt.plot(x0,y0,'bo',label='数据',color='black')
plt.plot(x0,[theta0+theta1*x for x in x0],label='批量梯度下降',color='red')
plt.plot(x0,[theta2+theta3*x for x in x0],label='随机梯度下降',color='blue')
plt.xlabel('x（面积）')
plt.ylabel('y（价格）')
plt.legend()
plt.show()
plt.scatter(range(count0),theta0_list,s=1)
plt.scatter(range(count0),theta1_list,s=1)
plt.xlabel('上方为theta0，下方为theta1')
plt.show()
plt.scatter(range(count1),theta2_list,s=3)
plt.scatter(range(count1),theta3_list,s=3)
plt.xlabel('上方为theta0，下方为theta1')
plt.show()

输出结果：

批量梯度下降最终得到theta0=2019.4181661097477，theta1=28.06952868382116
得到的回归函数是：y=2019.4181661097477+28.06952868382116x
随机梯度下降最终得到theta0=1685.432918583438，theta1=28.789137185954317
得到的回归函数是：y=1685.432918583438+28.789137185954317x
在这里插入图片描述
批量梯度下降时的 theta0 和 theta1 的变化：

随机梯度下降时的 theta0 和 theta1 的变化：

多元线性回归

1.定义数据

以下列数据为例
在这里插入图片描述
以上角标作为行索引，以下角标作为列索引

第二行可以写成：在这里插入图片描述
位于第二行第一列位置的数写成：
以上下角标来区分位置，以便于后期运算。

2.定义函数

设置一个回归方程：在这里插入图片描述
添加一个列向量

这样方程可以写为：在这里插入图片描述
既不会影响到方程的结果，而且使 $x$ 与 $\theta$ 的数量一致以便于矩阵计算。为了更方便表达，分别记为： $X=\begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \\ … \\ x_n \end{bmatrix},\theta=\begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \\ \theta_3 \\ … \\ \theta_n \end{bmatrix}$
这样方程就变为：
$h_\theta(x)=\begin{bmatrix} \theta_0 & \theta_1 & \theta_2 & \theta_3 & … & \theta_n\end{bmatrix}^T\begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \\ … \\ x_n \end{bmatrix}$
由回归方程推导出损失方程：
在这里插入图片描述

3，梯度下降

在这里插入图片描述
计算时要从 $\theta_0$ 一直计算到 $\theta_n$ 后，再从头由 $\theta_0$ 开始计算，以确保数据统一变化。在执行次数足够多的迭代后，我们就能取得达到要求的 $\theta_j$ 的值。

4鸢尾花数据集

Iris 鸢尾花数据集内包含 3 类，分别为山鸢尾（Iris-setosa）、变色鸢尾（Iris-versicolor）和维吉尼亚鸢尾（Iris-virginica），共 150 条记录，每类各 50 个数据，每条记录都有 4 项特征：花萼长度、花萼宽度、花瓣长度、花瓣宽度，可以通过这 4 个特征预测鸢尾花卉属于哪一品种。这是本文章所使用的鸢尾花数据集： sl：花萼长度；sw：花萼宽度；pl：花瓣长度；pw：花瓣宽度； type：类别：（Iris-setosa、Iris-versicolor、Iris-virginica 三类）在这里插入图片描述

import numpy as np
import pandas as pd
import random
import time
def MGD_train(X, y, alpha=0.0001, maxIter=1000, theta_old=None):
    '''
    MGD训练线性回归
    传入:
        X       :  已知数据  
        y       :  标签
        alpha   :  学习率
        maxIter :  总迭代次数 
    返回:
        theta : 权重参数
    '''
    # 初始化权重参数
    theta = np.ones(shape=(X.shape[1],))
    if not theta_old is None:
        # 假装是断点续训练
        theta = theta_old.copy()
    
    #axis=1 表示横轴，方向从左到右；axis=0 表示纵轴，方向从上到下
    for i in range(maxIter):
        # 预测
        y_pred = np.sum(X * theta, axis=1)
        # 全部数据得到的梯度
        gradient = np.average((y - y_pred).reshape(-1, 1) * X, axis=0)
        # 更新学习率
        theta += alpha * gradient
    return theta


def SGD_train(X, y, alpha=0.0001, maxIter=1000, theta_old=None):
    '''
    SGD训练线性回归
    传入:
        X       :  已知数据  
        y       :  标签
        alpha   :  学习率
        maxIter :  总迭代次数
        
    返回:
        theta : 权重参数
    '''
    # 初始化权重参数
    theta = np.ones(shape=(X.shape[1],))
    if not theta_old is None:
        # 假装是断点续训练
        theta = theta_old.copy()
    # 数据数量
    data_length = X.shape[0]
    for i in range(maxIter):
        # 随机选择一个数据
        index = np.random.randint(0, data_length)
        # 预测
        y_pred = np.sum(X[index, :] * theta)
        # 一条数据得到的梯度
        gradient = (y[index] - y_pred) * X[index, :]
        # 更新学习率
        theta += alpha * gradient
    return theta
    
    
def MBGD_train(X, y, alpha=0.0001, maxIter=1000, batch_size=10, theta_old=None):
    '''
    MBGD训练线性回归
    传入:
        X          :  已知数据  
        y          :  标签
        alpha      :  学习率
        maxIter    :  总迭代次数
        batch_size :  没一轮喂入的数据数
        
    返回:
        theta : 权重参数
    '''
    # 初始化权重参数
    theta = np.ones(shape=(X.shape[1],))
    
    if not theta_old is None:
        # 假装是断点续训练
        theta = theta_old.copy()
    
    # 所有数据的集合
    all_data = np.concatenate([X, y.reshape(-1, 1)], axis=1)
    for i in range(maxIter):
        # 从全部数据里选 batch_size 个 item
        X_batch_size = np.array(random.choices(all_data, k=batch_size))
        
        # 重新给 X, y 赋值
        X_new = X_batch_size[:, :-1]
        y_new = X_batch_size[:, -1]
        
        # 将数据喂入，更新 theta
        theta = MGD_train(X_new, y_new, alpha=0.0001, maxIter=1, theta_old=theta)
    return theta


def GD_predict(X, theta):
    '''
    用于预测的函数
    传入:
        X     : 数据
        theta : 权重
    返回:
        y_pred: 预测向量
    '''
    y_pred = np.sum(theta * X, axis=1)
    # 实数域空间 -> 离散三值空间，则需要四舍五入
    y_pred = (y_pred + 0.5).astype(int)
    return y_pred 


def calc_accuracy(y, y_pred):
    '''
    计算准确率
    传入:
        y        : 标签
        y_pred   : 预测值
    返回:
        accuracy : 准确率
    '''
    return np.average(y == y_pred)*100
    # 读取数据
iris_raw_data = pd.read_csv('iris.data', names  =['sepal length', 'sepal width', 'petal length', 'petal width', 'class'])

# 将三种类型映射成整数
Iris_dir = {'Iris-setosa': 1, 'Iris-versicolor': 2, 'Iris-virginica': 3}
iris_raw_data['class'] = iris_raw_data['class'].apply(lambda x:Iris_dir[x])


# 训练数据 X
iris_data = iris_raw_data.values[:, :-1]

# 标签 y
y = iris_raw_data.values[:, -1]

# 用 MGD 训练的参数
start = time.time()
theta_MGD = MGD_train(iris_data, y)
run_time = time.time() - start
y_pred_MGD = GD_predict(iris_data, theta_MGD)
print("MGD训练1000轮得到的准确率{:.2f}% 运行时间是{:.2f}s".format(calc_accuracy(y, y_pred_MGD), run_time))

# 用 SGD 训练的参数
start = time.time()
theta_SGD = SGD_train(iris_data, y)
run_time = time.time() - start
y_pred_SGD = GD_predict(iris_data, theta_SGD)
print("SGD训练1000轮得到的准确率{:.2f}% 运行时间是{:.2f}s".format(calc_accuracy(y, y_pred_SGD), run_time))

# 用 MBGD 训练的参数
start = time.time()
theta_MBGD = MBGD_train(iris_data, y)
run_time = time.time() - start
y_pred_MBGD = GD_predict(iris_data, theta_MBGD)
print("MBGD训练1000轮得到的准确率{:.2f}% 运行时间是{:.2f}s".format(calc_accuracy(y, y_pred_MBGD), run_time))