线性回归

线性回归:

如果只包括一个自变量和一个因变量,且二者的关系可用一条直线近似表示,这种回归分析称为一元线性回归分析;如果回归分析中包括两个或两个以上的自变量,且因变量和自变量之间是线性关系,则称为多元线性回归分析
线性回归是回归问题中的一种,线性回归假设目标值与特征之间线性相关,即满足一个多元一次方程

1 一元线性回归

1.1 什么是回归分析

回归分析是一种预测性的建模技术,它研究的是因变量(目标)和自变量(预测器)之间的关系。这种技术通常用于预测分析,时间序列模型以及发现变量之间的因果关系。通常使用曲线/线来拟合数据点,目标是使曲线到数据点的距离差异最小。

通过构建损失函数,来求解损失函数最小时的参数 w w w b b b。一元线性回归通常我们可以表达成如下公式:

y ^ \hat{y} y^为预测值,自变量 x x x和因变量 y y y是已知的,而我们想实现的是预测新增一个 x x x,其对应的 y y y是多少。因此,为了构建这个函数关系,目标是通过已知数据点,求解线性模型中 w w w b b b两个参数。

1.2 目标/损失函数

损失函数用来评价模型的预测值真实值不一样的程度,损失函数越好,通常模型的性能越好。不同的模型用的损失函数一般也不一样。在应用中,通常通过最小化损失函数求解和评估模型。

求解最佳参数,需要一个标准来对结果进行衡量,为此我们需要定量化一个目标函数式,使得计算机可以在求解过程中不断地优化。

针对任何模型求解问题,都是最终都是可以得到一组预测值 y ^ \hat{y} y^,对比已有的真实值 y y y ,数据行数为 n n n,可以将损失函数定义如下:

即预测值与真实值之间的平均的平方距离,统计中一般称其为MAE(mean square error)均方误差。把之前的函数式代入损失函数,并且将需要求解的参数w和b看做是函数L的自变量,可得:

现在的任务是求解最小化 L L L w w w b b b的值,

即核心目标优化式为:

求解方式:

1)最小二乘法(least square method)
不作为重点

2)梯度下降(gradient descent,GD)

梯度下降核心内容是对自变量进行不断的更新(针对w和b求偏导),使得目标函数不断逼近最小值的过程

2 梯度下降(GD)

2.1GD

img

2.2 核心公式

θ i = θ i − α ∂ ∂ θ i J ( θ ) θ_i=θ_i−α\frac{∂}{∂θ_i}J(θ) θi=θiαθiJ(θ)

其中 α α α步长,也叫学习率。步长(Learning rate)决定了在梯度下降迭代的过程中,每一步沿梯度负方向前进的长度。用上面下山的例子,步长就是在当前这一步所在位置沿着最陡峭最易下山的位置走的那一步的长度。
不同步长产生的效果不一样(如图所示):
因此,找到合适的步长十分重要。

2.3 梯度下降法的一般步骤

假设函数 y = f ( x 1 , x 2 , . . . , x n ) y = f (x_1,x_2,...,x_n) y=f(x1,x2,...,xn) 只有一个极小点。
初始给定参数为 X 0 = ( x 1 0 , x 2 0 , . . . , x n 0 , ) X_0 = (x_10,x_20,...,x_n0,) X0=(x10,x20,...,xn0,)。从这个点如何搜索才能找到原函数的极小值点?
方法:
①首先设定一个较小的正数 α α α ε ε ε
②求当前位置出处的各个偏导数:
f ′ ( x m 0 ) = ∂ y ∂ x m ( x m 0 ) , m = 1 , 2 , . . . , n f'(x_{m0})=\frac{∂y}{∂x_{m}}(x_{m0}),m=1,2,...,n f(xm0)=xmy(xm0),m=1,2,...,n
③修改当前函数的参数值,公式如下:
x m ′ = x m − α ∂ y ∂ x m ( x m 0 ) , m = 1 , 2 , . . . , n x'_{m}=x_{m}-α\frac{∂y}{∂x_{m}}(x_{m0}),m=1,2,...,n xm=xmαxmy(xm0),m=1,2,...,n
④如果参数变化量小于 ε ε ε,退出;否则返回第2步。

2.4 一元线性回归函数推导过程

设线性回归函数: y ^ = w x + b \hat{y}=wx+b y^=wx+b
构造损失函数(loss):
L ( w , b ) = 1 2 n ∑ i = 1 n ( w x i + b − y i ) 2 L(w,b)=\frac{1}{2n}\sum^n_{i=1}(wx_i+b-y_i)^2 L(w,b)=2n1i=1n(wxi+byi)2
思路:通过梯度下降法不断更新 w w w b b b,当损失函数的值特别小时,就得到了我们最终的函数模型。
过程:

step1.求导:
∂ ∂ w L ( w , b ) = 1 n ∑ i = 1 n ( ( w x i + b − y i ) ⋅ x i ) \frac{∂}{∂w}L(w,b)=\frac{1}{n}\sum^n_{i=1}((wx_i+b-y_i)·x_i) wL(w,b)=n1i=1n((wxi+byi)xi)

∂ ∂ b L ( w , b ) = 1 n ∑ i = 1 n ( w x i + b − y i ) \frac{∂}{∂b}L(w,b)=\frac{1}{n}\sum^n_{i=1}(wx_i+b-y_i) bL(w,b)=n1i=1n(wxi+byi)

step2.更新 θ 0 θ_{0} θ0 θ 1 θ_{1} θ1
w = w − α ⋅ 1 n ∑ i = 1 n ( ( w x i + b − y i ) ⋅ x i ) w=w-α·\frac{1}{n}\sum^n_{i=1}((wx_i+b-y_i)·x_i) w=wαn1i=1n((wxi+byi)xi)

b = b − α ⋅ 1 n ∑ i = 1 n ( w x i + b − y i ) b=b-α·\frac{1}{n}\sum^n_{i=1}(wx_i+b-y_i) b=bαn1i=1n(wxi+byi)

step3.代入损失函数,求损失函数的值,若得到的值小于 ε ε ε(一般为0.01或者0.001这样的小数),退出;否则,返回step1。

2.5 梯度下降存在的问题:

  • 当靠近极小值时收敛速度减慢;
  • 直线搜索时可能会产生一些问题;
  • 下降过程可能会出现“之字形”地下降。

2.6 三种梯度下降法(BGD,SGD,MBGD)

2.6.1 批量梯度下降法(Batch Gradient Descent)

批量梯度下降法每次都使用训练集中的所有样本更新参数。它得到的是一个全局最优解,但是每迭代一步,都要用到训练集所有的数据,如果m很大,那么迭代速度就会变得很慢。
优点:可以得出全局最优解。
缺点:样本数据集大时,训练速度慢。

2.6.2 随机梯度下降法(Stochastic Gradient Descent)

随机梯度下降法每次更新都从样本随机选择1组数据,因此随机梯度下降比批量梯度下降在计算量上会大大减少。SGD有一个缺点是,其噪音较BGD要多,使得SGD并不是每次迭代都向着整体最优化方向。而且SGD因为每次都是使用一个样本进行迭代,因此最终求得的最优解往往不是全局最优解,而只是局部最优解。但是大的整体的方向是向全局最优解的,最终的结果往往是在全局最优解附近。
优点:训练速度较快。
缺点:过程杂乱,准确度下降。

③小批量梯度下降法(Mini-batch Gradient Descent)

小批量梯度下降法对包含n个样本的数据集进行计算。综合了上述两种方法,既保证了训练速度快,又保证了准确度。

3 波士顿房价预测

3.1 代码

房屋价格与面积(数据在下面表格中):

序号面积价格
11506450
22007450
32508450
43009450
535011450
640015450
760018450

使用梯度下降求解线性回归(求 θ 0 θ_0 θ0 θ 1 θ_1 θ1
h θ ( x ) = θ 0 + θ 1 x h_θ(x)=θ_0+θ_1x hθ(x)=θ0+θ1x

#房屋价格与面积
#序号: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()

3.2 输出结果:

批量梯度下降最终得到theta0=2019.4181661097477,theta1=28.06952868382116
           得到的回归函数是:y=2019.4181661097477+28.06952868382116*x  
随机梯度下降最终得到theta0=1685.432918583438,theta1=28.789137185954317
           得到的回归函数是:y=1685.432918583438+28.789137185954317*x  

线性回归函数图像:

批量梯度下降时的theta0和theta1的变化:

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

4 多元线性回归及梯度下降

4.1 定义数据

以下面一组数据作为例子:

角标作为索引,以角标作为索引

第二行可以写成:
x ( 2 ) = [ 916 1201 5 … 33 ] x^{(2)}=\begin{bmatrix} 916 \\ 1201 \\ 5 \\ … \\ 33 \end{bmatrix} x(2)=9161201533
位于第二行第一列位置的数写成:
x 1 ( 2 ) = 916 x^{(2)}_1=916 x1(2)=916
以上下角标来区分位置,以便于后期运算。

4.2 多元线性回归 定义函数

设置一个回归方程:
h θ ( x ) = θ 0 + θ 1 x 1 + θ 2 x 2 + θ 3 x 3 + ⋯ + θ n x n h _\theta(x)=\theta_0+\theta_1x_1+\theta_2x_2+\theta_3x_3+\cdots+\theta_nx_n hθ(x)=θ0+θ1x1+θ2x2+θ3x3++θnxn
添加一个列向量
x 0 = [ 1 1 1 1 … 1 ] T x_0=\begin{bmatrix} 1 & 1 & 1 & 1 & … & 1\end{bmatrix}^T x0=[11111]T
这样方程可以写为:

h θ ( x ) = θ 0 x 0 + θ 1 x 1 + θ 2 x 2 + θ 3 x 3 + ⋯ + θ n x n h_\theta(x)=\theta_0x_0+\theta_1x_1+\theta_2x_2+\theta_3x_3+\cdots+\theta_nx_n hθ(x)=θ0x0+θ1x1+θ2x2+θ3x3++θnxn

既不会影响到方程的结果,而且使 x x x θ \theta θ的数量一致以便于矩阵计算。

为了更方便表达,分别记为:
X = [ x 0 x 1 x 2 x 3 … x n ] , θ = [ θ 0 θ 1 θ 2 θ 3 … θ n ] 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} X=x0x1x2x3xn,θ=θ0θ1θ2θ3θn
这样方程就变为:
h θ ( x ) = [ θ 0 θ 1 θ 2 θ 3 … θ n ] T [ x 0 x 1 x 2 x 3 … x n ] 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} hθ(x)=[θ0θ1θ2θ3θn]Tx0x1x2x3xn
由回归方程推导出损失方程:

4.3 梯度下降

计算时要从 θ 0 \theta_0 θ0一直计算到 θ n \theta_n θn后,再从头由 θ 0 \theta_0 θ0开始计算,以确保数据统一变化。

在执行次数足够多的迭代后,我们就能取得达到要求的 θ j \theta_j θj的值。

5 鸢尾花数据集

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

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-wLEijzbj-1618729225814)(https://jiang-hs.github.io/post-images/1594871083036.jpg)]

鸢尾花数据集下载:
https://archive.ics.uci.edu/ml/machine-learning-databases/iris/

下载这个iris.data即可
将其置于当前工作文件夹即可


先导入需要的库:

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'])

# 将三种类型映射成整数1.2.3
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))

运行结果:

MGD训练1000轮得到的准确率92.67% 运行时间是0.14s
SGD训练1000轮得到的准确率93.33% 运行时间是0.09s
MBGD训练1000轮得到的准确率92.67% 运行时间是0.14s

6 三种数据集

训练集、验证集、测试集。
通常比喻为学生的课本、作业及考试。
交叉验证
作用:防止过拟合。

思路为:利用训练集的数据进行参数评估,使用检验集为检验样本,计算并预测均方误差,比较各个模型的均方误差,选择误差最小的
评估模型是否学会了「某项技能」时,也需要用新的数据来评估,而不是用训练集里的数据来评估。这种「训练集」和「测试集」完全不同的验证方法就是交叉验证法。

7 回归模型评价指标

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

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值