机器学习算法的Python实现（一）：线性回归

在学习机器学习的过程中，结合数学推导和手写实现，可以加深对相关算法的认识。本部分教程将基于python实现机器学习的常用算法，来加强对算法的理解以及coding能力，仅供学习交流使用，请勿随意转载。

本篇内容从最基础的线性回归模型开始，全文分为三个部分：

• 数学推导
• python实现
• 线性回归优缺点分析

一、线性回归的数学推导

​ 在统计学中，线性回归（Linear Regression）是利用称为线性回归方程的最小平方函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析。这种函数是一个或多个称为回归系数的模型参数的线性组合。只有一个自变量的情况称为简单回归，大于一个自变量情况的叫做多元回归。

​ 在python中，可以通过以下8种方法实现线性回归：

1. Scipy.polyfit( ) or numpy.polyfit( )
2. Stats.linregress( )
3. Optimize.curve_fit( )
4. numpy.linalg.lstsq
5. Statsmodels.OLS ( )
6. 标准方程（简单的乘法求矩阵的逆）
7. 首先计算x的Moore-Penrose广义伪逆矩阵，然后与y取点积
8. sklearn.linear_model.LinearRegression( )

​ 线性模型就是对输入特征加权求和，再加上截距项，以此进行预测，如公式4-1所示。

公式1-1：线性回归模型预测

\[\hat{y}=\theta_0+\theta_1x_1+...+\theta_nx_n \]

​ 在此等式中：

• \(\hat{y}\) 是预测值
• \(n\) 是特征数量
• \(x_i\) 是第 \({i}\) 个特征值
• \(\theta_j\) 是第 \({j}\) 个模型参数（包括截距项 \(\theta_0\) 和特征权重 \(\theta_1,\theta_2,...,\theta_n\)）

可使用如公式1-2所示的向量化表示

公式1-2：线性回归模型预测（向量化表示）

\[\hat{y}=h_\theta(x)=\theta^{\mathrm{T}}\cdot x \]

​ 在此等式中：

• \(\theta\) 是模型的参数向量（列向量），其中包括截距项 \(\theta_0\) 和特征权重 \(\theta_1,\theta_2,...,\theta_n\)
• \(x\) 是模型的参数向量（列向量），包括 \(x_0, x_1,...,x_n\)，\(x_0\) 始终等于1
• \(\theta^{\mathrm{T}}\cdot x\) 是向量 \(\theta\) 和 \(x\) 的点积，其值等于 \(\theta_0x_0+\theta_1x_1+\theta_2x_2+...+\theta_nx_n\)
• \(h_\theta\) 是假设函数，使用模型参数 \(\theta\)

（一）损失函数——MSE

​ 如何训练回归模型使得模型参数拟合训练集并且具有较好的泛化性能呢。这里需要提到回归模型常用的性能指标——均方根误差（MSE），通过训练模型，获得最小化RMSE的 \(\theta\) 值。在实践中，为了求导方便，一般使用最小化均方误差（MSE）来达到训练效果。

公式1-3：线性回归模型的MSE函数

\[MSE=(X, h_\theta)=\frac{1}{m}\sum_{i=1}^{m}(\theta^{\mathrm{T}}x^{(i)}-y^{(i)})^2 \]

（二）求解MSE

​ 为了得到使成本函数最小的 \(\theta\) 值，可以使用最小二乘法获得标准方程求解或是梯度下降法进行迭代调整参数从而使成本函数最小化。先看一下这两种方法的对比：

	梯度下降法	正规方程法
学习速率 \(\alpha\)	需要设置	不需要
计算次数	需要多次迭代	不需要迭代
特征归一化	需要	不需要
时间复杂度	\(O(kn^2)\)	\(O(n^3)\)，需要计算 \(X^{\mathrm{T}}X\) 的逆
特征数量	即使n很大也可工作	如果n很大计算速度慢

1.标准方程

​ 要使 \(argmin(MSE)\)，可将其写成矩阵形式：

​ 假设 X ，是一个 m x (n+1) 的矩阵：每一行对应一个单独的训练样本

\[X=\left[ \begin{matrix} (x^{(1)})^{\mathrm{T}}\\ (x^{(2)})^{\mathrm{T}}\\ \ldots\\ (x^{(m)})^{\mathrm{T}}\\ \end{matrix} \right] \]

​ y 为一个 m 维的向量：包含所有训练集中的标签

\[y=\left[ \begin{matrix} y^{(1)}\\ y^{(2)}\\ \ldots\\ y^{(m)}\\ \end{matrix} \right] \]

​ 将MSE写成矩阵形式：\(MSE=\frac{1}{m}(X\theta-y)^{\mathrm{T}}(X\theta-y)\)，为了得到最小化，将MSE对 \(\theta\) 求导，结合矩阵求导的知识：

公式1-4：标准方程

\[\frac{\partial{MSE}}{\partial\theta}=\frac{\partial(\theta^{\mathrm{T}}X^{\mathrm{T}}X\theta-y^{\mathrm{T}}X\theta-\theta^{\mathrm{T}}X^{\mathrm{T}}y+y^{\mathrm{T}}y)}{\partial\theta}\\ \begin{align} \end{align} \]

​

​其中：


​第一项求导结果为 \(2X^{\mathrm{T}}X\theta\)


第二项求导结果为 \(-X^{\mathrm{T}}y\)
​第三项求导结果为 \(-X^{\mathrm{T}}y\)


第四项求导结果为0

​ 令其为0，则有：

\[\frac{\partial{MSE}}{\partial\theta} \begin{align} & =\frac{\partial(\theta^{\mathrm{T}}X^{\mathrm{T}}X\theta-y^{\mathrm{T}}X\theta-\theta^{\mathrm{T}}X^{\mathrm{T}}y+y^{\mathrm{T}}y)}{\partial\theta}\\ & = 2X^{\mathrm{T}}X\theta-2X^{\mathrm{T}}y=0\\ \end{align} \]

​ 如果 \(X^{\mathrm{T}}X\)可逆，则 \(\hat{\theta}=(X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}y\)，就是通过最小二乘法求解得到的使损失函数最小的 \(\theta\) 值。但是，当 \(X^{\mathrm{T}}X\)不可逆时，标准方程可能没有解。

2.最优化——梯度下降法

​ 梯度下降的中心思想是通过迭代调整参数，从而使损失函数最小化。根据在迭代过程中使用的数据集，分为批量梯度下降法、随机梯度下降法和小批量梯度下降法，我们这里使用的批量梯度下降（在应用梯度下降时，需要保证所有特征值大小比例都差不多，以尽可能减少收敛时间）。

​ 要实现梯度下降，就需要计算每个模型关于参数 \(\theta_j\) 的成本函数的梯度。单独计算每个模型的偏导数公式如下：

公式1-5：损失函数的偏导数

\[\frac{\partial{MSE(\theta)}}{\partial\theta_j}=\frac{2}{m}\sum_{i=1}^{m}(\theta^{\mathrm{T}}x^{(i)}-y^{(i)})x_j^{(i)} \]

​ 也可以使用公式1-6对所有成本函数的偏导数进行计算。

公式1-6：损失函数的梯度向量

\[\nabla{MSE}(\theta)= \left( \begin{matrix} \frac{\partial{MSE(\theta)}}{\partial\theta_0}\\ \frac{\partial{MSE(\theta)}}{\partial\theta_1}\\ \ldots\\ \frac{\partial{MSE(\theta)}}{\partial\theta_n} \end{matrix} \right) =\frac{2}{m}X^{\mathrm{T}}(X\theta-y) \]

公式1-7：梯度下降参数更新

\[\theta^{下一步}=\theta-\eta\nabla_{\theta}MSE(\theta) \]

梯度下降法的伪代码：

for 每一步 in 所有的训练批次:
使用整批训练数据获得梯度

​ 进行参数更新

二、Python实现

​ 在对线性回归的数学原理进行大致了解之后，我们开始进行算法的编写。我们首先构造数据集，然后根据梯度下降法的思路，我们需要在每个批次中对参数进行更新，这就需要我们在训练之初对参数进行初始化，然后根据损失函数获得损失函数的参数求导结果——即梯度，基于梯度对参数进行更新。最后我们使用交叉验证获得更加稳健的参数估计值。

​ 导入需要的相关包

import numpy as np
from sklearn.utils import shuffle
from sklearn.datasets import load_diabetes

（一）数据集构建

​ 本文基于sklearn自带的糖尿病数据集，进行简单说明。

def prepareData(self):
    """
    return:
    X-(442, 10)
    y-(442, 1)
    """
    # 调用sklearn的数据集接口获得相关数据
    data = load_diabetes().data
    target = load_diabetes().target
    # 打乱数据集
    X, y = shuffle(data, target, random_state=42)
    X, y = X.astype(np.float32), y.reshape((-1, 1))
    # 返回包含X与y的维度为（442，11）的数组
    data = np.c_[X, y]
    return data

（二）初始化参数

​ 根据输入数据的维度，返回维度分别为（dims，1）的参数w和（1，）的截距项。

def initializeParams(self, dims):
    """
    :param dims: X的维度10
    :return:
     w-(10, 1)
     b-(1,)
     """
    w = np.zeros((dims, 1))
    b = 0
    return w, b

（三）损失函数

​ 根据假设函数得到 \(\hat{y}\)，并计算损失函数，求出对应参数偏导数。

def linearLoss(self, X, y, w, b):
    """
    :param X: (442, 10)
    :param y: (10, 1)
    :param w: (10, 1)
    :param b: (1, )
    :return:
     y_hat-(442, 1)
     loss-(442, 1)
     dw-(10, 1)
     db-(1, )
     """
    num_train = X.shape[0]
    num_feature = 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 linearTrain(self, X, y, learning_rate, epochs):
    w, b = self.initializeParams(X.shape[1])
    losses = []
    for i in range(1, epochs+1):
        y_hat, loss, dw, db = self.linearLoss(X, y, w, b)
        losses.append(loss)
        # 梯度下降更新参数
        w += -learning_rate * dw
        b += -learning_rate * db

        # 打印迭代次数和损失
        if i % 10000 == 0:
            print(f'epoch:{i} loss: {loss}')

            # 保存参数
            params = {
                'w': w,
                'b': b
            }

            # 保存梯度
            grads = {
                'dw': dw,
                'db': db
            }

	return loss, params, grads

（五）交叉验证

​ 这里使用交叉验证获得测试集和验证集，返回的是一个生成器。

def linearCrossValidation(self, data, k, randomize=True):
	if randomize:
        data = list(data)
        shuffle(data)
    slices = [data[i::k] for i in range(k)]
    for i in range(k):
        validation = slices[i]
        train = [
            data
            for s in slices
            if s is not validation
            for data in s
        ]
        train = np.array(train)
        validation = np.array(validation)
        yield train, validation

（六）完整代码

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# author: GeekThomas
# time: 2021/3/21
import numpy as np
from sklearn.utils import shuffle
from sklearn.datasets import load_diabetes


class LrModel:
    def __init__(self):
        pass

    def prepareData(self):
        """
        return:
        X-(442, 10)
        y-(442, 1)
        """
        data = load_diabetes().data
        target = load_diabetes().target
        X, y = shuffle(data, target, random_state=42)
        X, y = X.astype(np.float32), y.reshape((-1, 1))
        data = np.c_[X, y]
        return data

    def initializeParams(self, dims):
        """
        :param dims: X的维度10
        :return:
        w-(10, 1)
        b-(1,)
        """
        w = np.zeros((dims, 1))
        b = 0
        return w, b

    def linearLoss(self, X, y, w, b):
        """
        :param X: (442, 10)
        :param y: (10, 1)
        :param w: (10, 1)
        :param b: (1, )
        :return:
        y_hat-(442, 1)
        loss-(442, 1)
        dw-(10, 1)
        db-(1, )
        """
        num_train = X.shape[0]
        num_feature = 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 linearTrain(self, X, y, learning_rate, epochs):
        w, b = self.initializeParams(X.shape[1])
        losses = []
        for i in range(1, epochs+1):
            y_hat, loss, dw, db = self.linearLoss(X, y, w, b)
            losses.append(loss)
            # 梯度下降更新参数
            w += -learning_rate * dw
            b += -learning_rate * db

            # 打印迭代次数和损失
            if i % 10000 == 0:
                print(f'epoch:{i} loss: {loss}')

            # 保存参数
            params = {
                'w': w,
                'b': b
            }

            # 保存梯度
            grads = {
                'dw': dw,
                'db': db
            }

        return loss, params, grads
	
    # 根据梯度下降法更新的参数，对模型进行预测，查看在测试集上的表现
    def predict(self, X, params):
        w, b = params['w'], params['b']
        y_pred = np.dot(X, w) + b
        return y_pred

    def linearCrossValidation(self, data, k, randomize=True):
        if randomize:
            data = list(data)
            shuffle(data)
        slices = [data[i::k] for i in range(k)]
        for i in range(k):
            validation = slices[i]
            train = [
                data
                for s in slices
                if s is not validation
                for data in s
            ]
            train = np.array(train)
            validation = np.array(validation)
            yield train, validation


if __name__ == '__main__':
    lr = LrModel()
    data = lr.prepareData()

    i = 1
    # 获得训练集和验证集
    for train, validation in lr.linearCrossValidation(data, 5):
        X_train, y_train = train[:, :10], train[:, -1].reshape((-1, 1))
        X_valid, y_valid = validation[:, :10], validation[:, -1].reshape((-1, 1))

        losses_5 = []
        loss, params, grads = lr.linearTrain(X_train, y_train, 0.001, 100000)
        losses_5.append(loss)
        score = np.mean(losses_5)
        print(f'{i} of 5 kold cross validation score is {score}')
        y_pred = lr.predict(X_valid, params)
        valid_score = np.sum(((y_pred - y_valid) ** 2)) / len(X_valid)
        print('valid score is', valid_score)
        i += 1

三、线性回归模型优缺点分析

（一）优点

1. 建模速度快，不需要很复杂的计算，在数据量大的情况下依然运行速度很快
2. 模型可解释性强。相对于很多机器学习的黑盒模型，线性回归可以根据系数给出每个变量的理解和解释

（二）缺点

​ 不能很好地拟合非线性数据。所以需要先判断变量之间是否是线性关系。

---

​ 以上就是本文的全部内容——基于Python实现简单线性回归。

参考资料：

1.机器学习：用正规方程法求解线性回归

2.矩阵求导

3.机器学习实验室：机器学习公式推导与代码实现