機器學習算法的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.機器學習實驗室：機器學習公式推導與代碼實現