07 模型之母：多元线性回归

2019-12-15 本文已影响0人 Japson

本文为饼干Japson原创，首发于公众号《数据科学家联盟》。更多大数据、机器学习、深度学习相关内容，敬请关注公众号。

0 前言

在线性回归的前3篇中，我们介绍了简单线性回归这种样本只有一个特征值的特殊形式，并且了解了一类机器学习的建模推导思想，即：

通过分析问题，确定问题的损失函数或者效用函数；

然后通过最优化损失函数或者效用函数，获得机器学习的模型。
然后我们推导并实现了最小二乘法，然后实现了简单线性回归。最后还以简单线性回归为例，学习了线性回归的评价指标：均方误差MSE、均方根误差RMSE、平均绝对MAE以及R方。

但是，在真实世界中，一个样本通常有很多（甚至成千上万）特征值的，这就是多元线性回归。本篇内容我们学习多元线性回归并实现。

1 多元线性回归

对于下面的样本数据集 $x^{(i)} = (X_1^{(i)},X_2^{(i)},...,X_n^{(i)})$ 对应的是一个向量，每一行是一个样本，每列对应一个特征。对应的结果可以用如下如下公式：
$y = \theta_0 + \theta_1x_1 + \theta_2x_2 + ... + \theta_nx_n$

简单线性回归，只计算前两项，但是在多元线性回归中就要学习到n+1个参数，就能求出多元线性回归预测值：
$\hat{y}_{(i)} = \theta_0 + \theta_1X_1^{(i)} + \theta_2X_2^{(i)} + ... + \theta_nX_n^{(i)}$

也就是：第一个特征与参数1相乘、第二个特征与参数2相乘，累加之后再加上截距。就能得到预测值。

求解思路也与简单线性回归非常一致，目标同样是：

已知训练数据样本 $x$ 、 $y$ ，找到 $\theta_0,\theta_1,\theta_2,...,\theta_n$ ，使 $\sum_{i=1}^{m} (y^{(i)}-\hat{y}^{(i)})^{2}$ 尽可能小.

其中 $\theta = (\theta_0,\theta_1,\theta_2,...,\theta_n)$ 是列向量列向量，而且我们注意到，可以虚构第0个特征X0，另其恒等于1，推导时结构更整齐，也更加方便：

$\hat{y}_{(i)} = \theta_0X_0^{(i)} + \theta_1X_1^{(i)} + \theta_2X_2^{(i)} + ... + \theta_nX_n^{(i)}$

这样我们就可以改写成向量点乘的形式：

15725963606023.jpg

此时，我们可以得出：

$\hat{y} = X_b · \theta$

因此我们可以把目标写成向量化的形式：

已知训练数据样本 $x$ 、 $y$ ，找到向量 $\theta$ ，使 $(y - X_b · \theta)^T(X_b · \theta)$ 尽可能小.

推导出可以得到多元线性回归的正规方程解：
$\theta = (X^T_bX_b)^{-1}X_b^{T}y$

当然了，具体的推导过程不需要了解的，不影响我们的使用，我们只要知道结果思想就行，结果也不用背下来，在网上搜一下就能找到。

但是这种朴素的计算方法，缺点是时间复杂度较高： $O(n^3)$ ，在特征比较多的时候，计算量很大。优点是不需要对数据进行归一化处理，原始数据进行计算参数，不存在量纲的问题（多选线性没必要做归一化处理）。

2 多元线性回归的实现

下面我们来使用python代码实现多元线性回归：

import numpy as np
from .metrics import r2_score

class LinearRegression:

    def __init__(self):
        """初始化Linear Regression模型"""
        self.coef_ = None    # 系数（theta0~1 向量）
        self.interception_ = None   # 截距（theta0 数）
        self._theta = None  # 整体计算出的向量theta

    def fit_normal(self, X_train, y_train):
        """根据训练数据X_train，y_train训练Linear Regression模型"""
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"
        # 正规化方程求解
        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]
        return self

    def predict(self, X_predict):
        """给定待预测的数据集X_predict，返回表示X_predict的结果向量"""
        assert self.interception_ is not None and self.coef_ is not None, \
            "must fit before predict"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"
        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        y_predict = X_b.dot(self._theta)
        return y_predict

    def score(self, X_test, y_test):
        """很倔测试机X_test和y_test确定当前模型的准确率"""
        y_predict = self.predict(self, X_test)
        return r2_score(y_test, y_predict)
    

    def __repr__(self):
        return "LinearRegression()"

其实在代码中，思想很简单，就是使用公式即可。其中有一些知识点：

1、np.hstack(tup)：参数tup可以是元组，列表，或者numpy数组，返回结果为numpy的数组。按列顺序把数组给堆叠起来（加一个新列）。

2、np.ones()：返回一个全1的n维数组，有三个参数：shape（用来指定返回数组的大小）、dtype（数组元素的类型）、order（是否以内存中的C或Fortran连续（行或列）顺序存储多维数据）。后两个参数都是可选的，一般只需设定第一个参数。（类似的还有np.zeros()返回一个全0数组）

3、numpy.linalg模块包含线性代数的函数。使用这个模块，可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等。inv函数计算逆矩阵

4、T：array的方法，对矩阵进行转置。

5、dot：点乘

3 调用

下面我们可以在jupyter notebook中调用我们的算法：

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

boston = datasets.load_boston()

X = boston.data
y = boston.target

X = X[y<50.0]
y = y[y<50.0]

X.shape
输出：(490, 13)

y.shape
输出：(490, )

from myAlgorithm.model_selection import train_test_split
from myAlgorithm.LinearRegression import LinearRegression

X_train, X_test, y_train, y_test = train_test_split(X, y, seed = 666)

reg = LinearRegression()
reg.fit_normal(X_train, y_train)

reg.coef_
输出：
array([-1.18919477e-01,  3.63991462e-02, -3.56494193e-02,  5.66737830e-02,
       -1.16195486e+01,  3.42022185e+00, -2.31470282e-02, -1.19509560e+00,
        2.59339091e-01, -1.40112724e-02, -8.36521175e-01,  7.92283639e-03,
       -3.81966137e-01])

reg.interception_
输出：
34.16143549622471

reg.score(X_test, y_test)
输出：
0.81298026026584658

我们看到，reg.coef_这一项的结果是13个系数，这13个系数有正有负。正负代表着该系数所乘的特征与预测目标是正相关还是负相关。正相关，特征越大房价越高；负相关，特征越大，房价越低。而系数绝对值的大小决定了影响程度。

下面我们对所有的系数按照数值由小到大进行排序：

np.argsort(reg.coef_)
输出：
array([ 4,  7, 10, 12,  0,  2,  6,  9, 11,  1,  3,  8,  5])

将这个返回结果作为索引，返回排序后索引所对应的特征名：

boston.feature_names[np.argsort(reg.coef_)]
输出：
array(['NOX', 'DIS', 'PTRATIO', 'LSTAT', 'CRIM', 'INDUS', 'AGE', 'TAX',
       'B', 'ZN', 'CHAS', 'RAD', 'RM'], dtype='<U7')

这也说明了线性回归算法，具有可解释性。

4 总结

在本篇内容中，首先学习了多元线性回归的推导过程，然后使用代码实现，最后进行了调用。

到此为止，线性回归模型就介绍完了。线性回归模型有着比较清晰的数据推导过程，也是其他复杂模型的基础。线性回归算法是典型的参数学习。虽然线性回归只能解决回归问题，但是却是很多分类问题，如逻辑回归的基础。并且线性回归算法是假设数据是有一定的线性关系的，且线性关系越强，效果越好。

在第一节中得到的多元线性回归的正规方程解，看上去很简单，但是时间复杂度高。其实除了使用正规方程解以外，还可以使用大名鼎鼎的梯度下降法。梯度下降法不仅可以解决线性问题，更是解决机器学习的最优模型的通用算法。So，下面就是梯度下降的学习啦。

公众号二维码.jpeg