Press "Enter" to skip to content

手撸机器学习算法 – 线性回归

本站内容均来自兴趣收集,如不慎侵害的您的相关权益,请留言告知,我们将尽快删除.谢谢.

系列文章目录:

感知机
线性回归

如果说感知机是最最最简单的分类算法,那幺线性回归就是最最最简单的回归算法,所以这一篇我们就一起来快活的用两种姿势手撸线性回归吧;

 

算法介绍

 

线性回归通过超平面拟合数据点,经验误差一般使用MSE(均平方误差),优化方法为 最小二乘法 ,算法如下:

 

 

    1. 假设输入数据为X,输出为Y,为了简单起见,这里的数据点为一维数据(更好可视化,处理方式没区别);

 

    1. MSE公式为: \(\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\) ;

 

    1. 最小二乘法:最小指的是目标是 min ,二乘指的就是MSE中 误差的二次方 ,公式为: \(min\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\) ;

 

    1. 由于目标是查找拟合最好的超平面,因此依然定义变量 w 和 b ;

 

    1. 对于w和b的求解有 两种 方式:

 

      1. 列出最小化的公式,利用优化求解器求解:

 

        1. 基于已知的X、Y,未知的w、b构建MSE公式;

 

        1. 定义最小化MSE的目标函数;

 

        1. 利用求解器直接求解上述函数得到新的w和b;

 

对经验误差函数求偏导并令其为0推导出 w 和 b 的解析解:

 

        1. 基于最小化MSE的优化问题可以直接推导出w和b的计算方法;

 

        1. 基于推导出的计算方法直接计算求解;

 

 

 

 

利用求解器求解

 

利用求解器求解可以看作就是个列公式的过程,把已知的数据X和Y,未知的变量w和b定义好,构建出MSE的公式,然后丢到求解器直接对w和b求偏导即可,相对来说代码繁琐,但是过程更简单,没有任何数学推导;

 

代码实现

 

初始化数据集

 

X = np.array([1.51, 1.64, 1.6, 1.73, 1.82, 1.87])
y = np.array([1.63, 1.7, 1.71, 1.72, 1.76, 1.86])

 

定义变量符号

 

所谓变量指的就是那些需要求解的部分,次数就是超平面的w和b;

 

w,b = symbols('w b',real=True)

 

定义经验误差函数MSE

 

RDh = 0
for xi,yi in zip(X,y):
    RDh += (yi - (w*xi+b))**2
RDh = RDh / len(X)

 

定义求解函数

 

此处就是对w和b求偏导;

 

eRDHw = diff(RDh,w)
eRDHb = diff(RDh,b)

 

求解w和b

 

ans = solve((eRDHw,eRDHb),(w,b))
w,b = ans[w],ans[b]

 

运行结果

 

 

完整代码

 

from sympy import symbols, diff, solve
import numpy as np
import matplotlib.pyplot as plt
'''
线性回归拟合wx+b直线;
最小二乘法指的是优化求解过程是通过对经验误差(此处是均平方误差)求偏导并令其为0以解的w和b;
'''
# 数据集 D X为父亲身高,Y为儿子身高
X = np.array([1.51, 1.64, 1.6, 1.73, 1.82, 1.87])
y = np.array([1.63, 1.7, 1.71, 1.72, 1.76, 1.86])
# 构造符号
w,b = symbols('w b',real=True)
# 定义经验误差计算公式:(1/N)*sum(yi-(w*xi+b))^2)
RDh = 0
for xi,yi in zip(X,y):
    RDh += (yi - (w*xi+b))**2
RDh = RDh / len(X)
# 对w和b求偏导:求偏导的结果是得到两个结果为0的方程式
eRDHw = diff(RDh,w)
eRDHb = diff(RDh,b)
# 求解联立方程组
ans = solve((eRDHw,eRDHb),(w,b))
w,b = ans[w],ans[b]
print('使得经验误差RDh取得最小值的参数为:'+str(ans))
plt.scatter(X,y)
x_range = [min(X)-0.1,max(X)+0.1]
y_range = [w*x_range[0]+b,w*x_range[1]+b]
plt.plot(x_range,y_range)
plt.show()

 

推导公式求解

 

与利用优化器求解的区别在于针对 \(min\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\) 对 \(w\) 和 \(b\) 求偏导并令其为0,并推导出 w 和 b 的计算公式是自己推导的,还是由优化器完成的,事实上如果自己推导,那幺最终代码实现上会非常简单(推导过程不会出现在代码中);

 

w和b的求解公式推导

 

首先,我们的优化目标为:

 

\[min \frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2 \]

 

去除公式中无关的常量部分:

 

\[min \sum_{i=1}^{N}(w*x_i+b-y_i)^2 \]

 

由于一般w是向量,而b为标量,因此通常会将w和b组成[w b],x变为[x 1]来统一处理w和b,调整后如下:

 

\[\sum_{i=1}^{N}(wx_i^T-y_i)^2 \]

 

上式把平方拆开有:

 

\[\sum_{i=1}^{N}(ww^Tx_ix_i^T-2wx_i^Ty_i+y_i^2) \]

 

对于w(注意此时w为原w和b的组合)求偏导过程如下:

 

\[\begin{equation*} \begin{split} \frac{\partial \sum_{i=1}^{N}(ww^Tx_ix_i^T-2wx_i^Ty_i+y_i^2)}{\partial w} &= 2w^Tx_ix_i^T-2x_i^Ty_i \\ &= 0 \\ \end{split} \end{equation*} \]

 

上式变形后有:

 

\[2w^Tx_ix_i^T – 2x_i^Ty_i = 0 \\ w^Tx_ix_i^T = x_i^Ty_i \\ w^T = (x_ix_i^T)^{-1}x_i^Ty_i \]

 

由于此处的w其实是w和b的组合,因此通过这一次推导就得到了w和b两个求解方法;

 

代码实现

 

构造数据集

 

X = np.array([1.51,1.64,1.6,1.73,1.82,1.87]).reshape(-1,1)
y = np.array([1.63,1.7,1.71,1.72,1.76,1.86])

 

为X增加元素全为1的一列用于和b的计算

 

ones = np.ones(X.shape[0]).reshape(-1,1)
X = np.hstack((ones,X))

 

通过求解公式求解w和b

 

w = np.linalg.inv(self.X.T @ self.X) @ self.X.T @ self.y
w,b = w[1:],w[0]

 

运行结果

 

 

完整代码

 

完整代码对于求解部分使用的是伪逆而不是逆,原因在于求解公式中正好构造了伪逆,而伪逆适用性强国求逆,因此使用伪逆代替逆;

 

import numpy as np
import matplotlib.pyplot as plt
rnd = np.random.RandomState(3)  # 为了演示,采用固定的随机
'''
单变量线性回归最小二乘法的矩阵实现:矩阵实现的优势在于numpy本身支持伪逆;
其实就是对于误差平方和的矩阵形式对于W求导并令其为0,得到w_hat = (X^T*X)^-1*X^T*Y,其中(X^T*X)^-1*X^T称为伪逆(pseudo inverse,即函数pinv)
因此可以省略中间大量的构造经验误差、解偏导方程组等步骤;
'''
class LinearRegression(object):
    def __init__(self,X,y):
        ones = np.ones(X.shape[0]).reshape(-1,1) # 1用于计算b
        self.X = np.hstack((ones,X))
        self.y = y
    def train(self):
        # 注意,虽然一般情况下下面二者是等价的,但是在矩阵无法求逆或某些其他情况下时,二者并不相等
        # 相对而言伪逆定义更加宽泛,用处更广,因此可以的情况下建议使用伪逆
        # self.w = np.linalg.inv(self.X.T @ self.X) @ self.X.T @ self.y
        self.w = np.linalg.pinv(self.X) @ self.y
        self.w = self.w.reshape(-1)
        self.w,self.b = self.w[1:],self.w[0]
        return self.w,self.b
    def predict(self,x):
        return self.w.dot(x)+self.b
    def get(self):
        return self.X,self.y,self.w,self.b
if __name__ == '__main__':
    X0 = np.array([1.51,1.64,1.6,1.73,1.82,1.87]).reshape(-1,1)
    y = np.array([1.63,1.7,1.71,1.72,1.76,1.86])
    model = LinearRegression(X=X0,y=y)
    w,b = model.train()
    print(f'最小二乘法的矩阵方式结果为:w={w} b={b}')
    print(model.predict(np.array([X0[0]])))
    
    plt.scatter(X0,y)
    plt.plot([min(X0),max(X0)],[model.predict(min(X0)),model.predict(max(X0))])
    plt.show()

 

最后

 

对于算法的学习,一定的数学知识是必要的,对于公式的推导可以让我们对于算法内部运行逻辑有更深的了解;

Be First to Comment

发表评论

您的电子邮箱地址不会被公开。