线性回归（公式推导，Numpy、sklearn实现）

本文章使用的数据集以及源码Github-stellaris123

线性回归

​ 线性回归是线性模型的一种典型方法。产品销量预测、岗位薪资预测，都可以用先线性回归来拟合模型。从某种程度上来说，回归分析不再局限于线性回归这一具体模型和算法，更包含了广泛的由自变量到因变量的机器学习建模思想。

原理

​ 给定一组由输入 $x$ 和输出 $y$ 构成的数据集 $D=(x_1,y_1),\dots ,(x_k,y_k)$，其中 $x_i$ 为参数集合。线性回归就是通过不断训练从而得到一个线性模型去尽可能的根据输入 $x$ 拟合出 $y$ 值。

​ 以 $x$ 为影响因素，$y$ 为输出结果，构建回归式：$y=w{x}_i+b$，其中 $w$ 为模型矩阵。

​ 线性回归模型的关键问题是确定 $w$ 和 $b$ 的值，使得拟合输出 $y$ 与真实值 $y_i$ 尽可能接近。在回归任务中，我们通常使用均方误差来度量预测与标签值之间的损失，所以回归任务的优化目标就是使得拟合输出和真实输出之间的均方误差最小化。
$$\begin{array}{rl}f(w^{\ast },b^{\ast })& =argmin\sum _{i=1}^k(y-y_i{)}^2\\ & =argmin\sum _{i=1}^k(w{x}_i+b-y_i{)}^2\end{array}\text{\hspace{1em}(1)}$$
​ 为求的 $w$ 和 $b$ 的最小化参数 $w^{\ast }$ 和 $b^{\ast }$，可从式 $(1)$，分别对 $w$ 和 $b$ 求一阶导为零。

对 $w$ 求导
$$\begin{array}{rl}\frac{\partial f(w,b)}{\partial w}& =\frac{\partial }{\partial w}[\sum _{i=1}^k(w{x}_i+b-y_i{)}^2]\\ & =\sum _{i=1}^k\frac{\partial }{\partial w}[(y_i-w{x}_i-b{)}^2]\\ & =\sum _{i=1}^k[2\cdot (y_i-w{x}_i-b)\cdot (-x_i)]\\ & =\sum _{i=1}^k[2\cdot (w{x}_i^2-y_i{x}_i+b{x}_i)]\\ & =2\cdot (w\sum _{i=1}^k{x}_i^2-\sum _{i=1}^k{y}_i{x}_i+\sum _{i=1}^{k}b{x}_i)\end{array}\text{\hspace{1em}(2)}$$
对 $b$ 求导
$$\begin{array}{rl}\frac{\partial f(w,b)}{\partial b}& =\frac{\partial }{\partial b}[\sum _{i=1}^k(w{x}_i+b-y_i{)}^2]\\ & =\sum _{i=1}^k\frac{\partial }{\partial b}[(y_i-w{x}_i-b{)}^2]\\ & =\sum _{i=1}^k[2\cdot (y_i-w{x}_i-b)\cdot (-1)]\\ & =\sum _{i=1}^k[2\cdot (b-y_i+w{x}_i)]\\ & =2\cdot (\sum _{i=1}^{k}b-\sum _{i=1}^k{y}_i+\sum _{i=1}^{k}w{x}_i)\\ & =2\cdot (kb-\sum _{i=1}^k{y}_i+\sum _{i=1}^{k}w{x}_i)\end{array}\text{\hspace{1em}(3)}$$
令式 $(2)$ 为零
$$\begin{array}{rl}\frac{\partial f(w,b)}{\partial w}& =2\cdot (w\sum _{i=1}^k{x}_i^2-\sum _{i=1}^k{y}_i{x}_i+\sum _{i=1}^{k}b{x}_i)=0\\ & \iff w\sum _{i=1}^k{x}_i^2=\sum _{i=1}^k{y}_i{x}_i-\sum _{i=1}^{k}b{x}_i\end{array}\text{\hspace{1em}(4)}$$
令式 $(3)$ 为零
$$\begin{array}{rl}\frac{\partial f(w,b)}{\partial b}& =2\cdot (kb-\sum _{i=1}^k{y}_i+\sum _{i=1}^{k}w{x}_i)=0\\ & \iff b=\frac{1}{k}\sum _{i=1}^k(y_i-w{x}_i)\end{array}\text{\hspace{1em}(5)}$$
又因为：
$$\begin{gathered} \frac{1}{k}\sum _{i=1}^k{y}_i=\bar{y} \\ \frac{1}{k}\sum _{i=1}^k{x}_i=\bar{x}\text{\hspace{1em}(6)} \end{gathered}$$
所以：
$$b=\bar{y}-w\bar{x}\text{\hspace{1em}(7)}$$
将 $(7)$ 代入 $(4)$ 得
$$\begin{array}{rl}& w\sum _{i=1}^k{x}_i^2=\sum _{i=1}^k{y}_i{x}_i-\sum _{i=1}^{k}b{x}_i\\ & w\sum _{i=1}^k{x}_i^2=\sum _{i=1}^k{y}_i{x}_i-\sum _{i=1}^k(\bar{y}-w\bar{x})x_i\\ & w\sum _{i=1}^k{x}_i^2=\sum _{i=1}^k{y}_i{x}_i-\bar{y}\sum _{i=1}^k{x}_i+w\bar{x}\sum _{i=1}^k{x}_i\\ & w(\sum _{i=1}^k{x}_i^2-\bar{x}\sum _{i=1}^k{x}_i)=\sum _{i=1}^k{y}_i{x}_i-\bar{y}\sum _{i=1}^k{x}_i\\ & w=\frac{{\sum }_{i=1}^k{y}_i{x}_i-\bar{y}{\sum }_{i=1}^k{x}_i}{{\sum }_{i=1}^k{x}_i^2-\bar{x}{\sum }_{i=1}^k{x}_i}\end{array}\text{\hspace{1em}(8)}$$
又因为：
$$\begin{gathered} \bar{y}\sum _{i=1}^k{x}_i=\frac{1}{k}\sum _{i=1}^k{y}_i\sum _{i=1}^k{x}_i=\bar{x}\sum _{i=1}^k{y}_i \\ \bar{x}\sum _{i=1}^k{x}_i=\frac{1}{k}\sum _{i=1}^k{x}_i\sum _{i=1}^k{x}_i=\frac{1}{k}(\sum _{i=1}^k{x}_i{)}^2\text{\hspace{1em}(9)} \end{gathered}$$
所以：
$$w=\frac{{\sum }_{i=1}^k{y}_i(x_i-\bar{x})}{{\sum }_{i=1}^k{x}_i^2-\frac{1}{k}({\sum }_{i=1}^k{y}_i{)}^2}\text{\hspace{1em}(10)}$$

其他转化：我们可能会见到类似式 $(13)$ 的公式，这里说一下是怎么推导的

将 $(9)$ 代回 $(10)$
$$\begin{array}{rl}w& =\frac{{\sum }_{i=1}^k{y}_i(x_i-\bar{x})}{{\sum }_{i=1}^k{x}_i^2-\bar{x}{\sum }_{i=1}^k{x}_i}\\ & =\frac{{\sum }_{i=1}^k(y_i{x}_i-y_i\bar{x})}{{\sum }_{i=1}^k(x_i^2-x_i\bar{x})}\end{array}\text{\hspace{1em}(11)}$$
又因为：
$$\begin{gathered} \bar{y}\sum _{i=1}^k{x}_i=\bar{x}\sum _{i=1}^k{y}_i=\sum _{i=1}^k\bar{y}{x}_i=\sum _{i=1}^k\bar{x}{y}_i=k\bar{x}\bar{y}=\sum _{i=1}^k\bar{x}\bar{y} \\ \bar{x}\sum _{i=1}^k{x}_i=\sum _{i=1}^k\bar{x}{x}_i=k\bar{x}\bar{x}=\sum _{i=1}^k{\bar{x}}^2\text{\hspace{1em}(12)} \end{gathered}$$
所以：（这个公式也是常见的手推最小二乘法公式了）
$$\begin{array}{rl}w& =\frac{{\sum }_{i=1}^k(y_i{x}_i-y_i\bar{x}-x_i\bar{y}+\bar{x}\bar{y})}{{\sum }_{i=1}^k(x_i^2-x_i\bar{x}-x_i\bar{x}+{\bar{x}}^2)}\\ & =\frac{{\sum }_{i=1}^k(x_i-\bar{x})(y_i-\bar{y})}{{\sum }_{i=1}^k(x_i-\bar{x}{)}^2}\end{array}\text{\hspace{1em}(13)}$$

然后 $b$ 就可以用式 $(7)$ 来求了。

实现转换

多元线性回归

​ 多元问题，就是输入多个变量，也就是一个向量变量，这个时候对应的 $w$ 也是向量。这个是时候需要将上述的式子矩阵化。

如果想要用 Python 来实现的话，$(10)$ 中的求和运算只能用循环来实现。但是如果能将 $(10)$ 向量化，也就是转换成矩阵运算的话，就可以用 $Numpy$ 这种专门加速矩阵运算的类库来进行编写。

令：【$x,x_d,y,y_d$ 均为 $k$ 行 $1$ 列的列向量，$x_d,y_d$ 分别为去均值后的 $x,y$】
$$\begin{gathered} x=(x_1;x_2;\cdots ;x_k),x_d=(x_1-\bar{x};x_2-\bar{x};\cdots ;x_k-\bar{x}) \\ y=(y_1;y_2;\cdots ;y_k),y_d=(y_1-\bar{y};y_2-\bar{y};\cdots ;y_k-\bar{y})\text{\hspace{1em}(14)} \end{gathered}$$
将 $(14)$ 代入 $(13)$
$$w=\frac{x_d^T{y}_d}{x_d^T{x}_d}$$

Nnmpy实现

定义回归模型

import numpy as np
def linear_loss(X, y, w, b):
    """
    输入:
        X->变量矩阵
        y->标签向量
        w->权重矩阵
        b->偏置
    输出:
        y_hat->
    """
    # 训练样本量
    num_train = X.shape[0]
    # 训练特征量
    num_feawture = X.shape[1]
    # 线性回归预测值
    y_hat = np.dot(X, w) + b;
    # 计算均方损失
    loss = np.sum((y_hat - y)**2)/num_train
    # 基于均方损失对权重系数的一阶梯度
    dw = np.dot(X.T, (y_hat - y))/num_train
    # 基于均方损失对偏置的一阶梯度
    db = np.sum((y_hat - y))/num_train
    return y_hat, loss, dw, db

初始化模型参数

def init_params(dims):
    # 初始化权重系数为零向量
    w = np.zeros((dims,1))
    # 初始化偏置参数为零
    b = 0
    return w, b

迭代训练

def linear_train(X, y, learning_rate = 0.01, train_steps = 10000):
    """
    输入:
        X->变量矩阵
        y->标签向量
        learning_rate->学习率（默认0.01）
        train_steps->训练次数（默认10000）
    输出:
        loss_his->当前迭代的损失
        parmes->当前迭代步优化后的参数
        grads->当前迭代步的梯度
    """
    # 记录训练损失
    loss_his = []
    # 初始化模型参数
    w, b = init_params(X.shape[1])
    # 迭代训练
    for i in range(1, train_steps):
        # 计算当前迭代的预测值
        y_hat, loss, dw, db = linear_loss(X, y, w, b)
        # 基于梯度下降法的参更新
        w += -learning_rate * dw
        b += -learning_rate * db
        # 记录当前迭代的损失
        loss_his.append(loss)
        # 每10000次迭代打印当前损失信息
        if (i+1)%10000==0:
            print("epoch %d loss %f"%(i+1, loss))
        # 将当前迭代步优化后的参数保存到字典中
        params={
            'w':w,
            'b':b
        }
        # 将当前迭代步的梯度保存到字典中
        grads={
            'dw':dw,
            'db':db
        }
    return loss_his, params, grads

预测函数

def predict(X, params):
    """
    输入:
        X->测试数据集
        params->模型训练参数
    输出:
        y_pred->模型预测结果
    """
    # 获取模型参数
    w = params['w']
    b = params['b']
    # 预测
    y_pred = np.dot(X, w) + b
    return y_pred

R2系数函数

def r2_score(y_test, y_pred):
    """
    输入:
        y_test->测试集标签值
        y_pred->测试集预测值
    输出:
        r2->R2系数
    """
    # 测试标签均值
    y_avg = np.mean(y_test)
    # 总离差平方和
    ss_tot = np.sum((y_test - y_avg)**2)
    # 残差平方和
    ss_res = np.sum((y_test - y_pred)**2)
    # R2计算
    r2 = 1 - (ss_res/ss_tot)
    return r2

尾注

参考

[1]鲁伟.机器学习-公式推导与代码实现[M].北京：人民邮电出版社，2022.

[2]周志华.机器学习[M].北京：清华大学出版社，2016.

[3]谢文睿，秦川.机器学习公式详解[M].北京：人民邮电出版社，2022.

[4]李航.统计学习方法[M].北京：清华大学出版社，2019.