机器学习之线性回归法

摘要：在统计学中，线性回归是利用称为线性回归方程的最小二乘函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析维基百科。对当前示例作图表示衡量线性回归法的指标误差一个训练后的模型通常都会使用测试数据集测试该模型的准确性。

在统计学中，线性回归（Linear regression）是利用称为线性回归方程的最小二乘函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析^维基百科。

当只有一个自变量的时候，成为简单线性回归。

为了得到一个简单线性回归模型，假设存在以房屋面积为特征，以价格为样本输出，包含四个样本的样本集，如图：

寻找一条直线，最大程度上拟合样本特征与样本输出之间的关系。

假设最佳拟合的直线方程为：$y = ax + b$，则对于样本特征 $x$ 的每一个取值 $x^{(i)}$ 的预测值为：$hat{y}^{(i)} = ax^{(i)} +b$。而我们希望的就是真值 $y^{(i)}$ 和预测值 $hat{y}^{(i)}$ 之间的差距尽量小。

可以用 $|y^{(i)} - hat{y}^{(i)}|$ 表示两者之间的差距，对于所有的样本，使用求和公式求和处理：

$$ sum_{i=1}^{m}|y^{(i)} - hat{y}^{(i)}| $$

但是这个公式有一个问题，不容易求导，为了解决这个问题，可先对 $|y^{(i)} - hat{y}^{(i)}|$ 进行平方，如此最后的公式就变成了：

$$ sum_{i=1}^{m}(y^{(i)} - hat{y}^{(i)})^2 $$

最后，替换掉 $hat{y}^{(i)}$ ，即为：

$$ sum_{i=1}^{m}(y^{(i)} - ax^{(i)}-b)^2 $$

因此，找到的一个简单线性回归模型就是找到合适的 a 和 b，使得该函数的值尽可能的小，该函数也称为损失函数（loss function）。

找到合适的 a 和 b，使得 $sum_{i=1}^{m}(y^{(i)} - ax^{(i)}-b)^2$ 的值尽可能的小，这样的方法称为最小二乘法。

如何求 a 和 b 呢？令该函数为 $J(a, b)$，分别使对 a 和 b 求导的结果为0。

对 b 求导：$frac{partial{J(a, b)}}{partial{b}} = sum_{i=1}^{m}2(y^{(i)} - ax^{(i)}-b)(-1) = 0$，得：

$$ b = overline{y}-aoverline{x} $$

对 a 求导：$frac{partial{J(a, b)}}{partial{a}} = sum_{i=1}^{m}2(y^{(i)} - ax^{(i)}-b)(-x^{(i)}) = 0$，得：

$$ a = frac{sum_{i=1}^{m}(x^{(i)}-overline{x})(y^{(i)}-overline{y})}{sum_{i=1}^{m}(x^{(i)}-overline{x})^2} $$

注：这里略去了公式的推导过程。

有了数学的帮助，实现简单线性回归就比较方便了。

首先声明一个样本集：

import numpy as np

x = np.array([1., 2., 3., 4., 5.])
y = np.array([1., 3., 2., 3., 5.])

公式中用到了 x 和 y 的均值：

x_mean = np.mean(x)
y_mean = np.mean(y)

求 a 和 b 的值有两种方法。第一种是使用 for 循环：

# 分子
num = 0.0

# 分母
d = 0.0

for x_i, y_i in zip(x, y):
    num += (x_i - x_mean) * (y_i - y_mean)
    d += (x_i - x_mean) ** 2
    
a = num / d
b = y_mean - a * x_mean

第二种是使用矩阵乘：

num = (x - x_mean).dot(y - y_mean)
d = (x - x_mean).dot(x - x_mean)

a = num / d
b = y_mean - a * x_mean 

注：使用矩阵乘效率更高。

求出了 a 和 b，简单线性模型就有了：$hat{y} = a*x + b$。对当前示例作图表示：​

一个训练后的模型通常都会使用测试数据集测试该模型的准确性。对于简单线性归回模型当然可以使用 $sum_{i=1}^{m}(y_{test}^{(i)} - hat{y}_{test}^{(i)})^2$ 来衡量，但是它的取值和测试样本个数 m 存在联系，改进方法很简单，只需除以 m 即可，即均方误差（Mean Squared Error）：

$$ MSE：frac{1}{m}sum_{i=1}^{m}(y_{test}^{(i)} - hat{y}_{test}^{(i)})^2 $$

np.sum((y_predict - y_true) ** 2) / len(y_true)

值得一提的是 MSE 的量纲是样本单位的平方，有时在某些情况下这种平方并不是很好，为了消除量纲的不同，会对 MSE 进行开方操作，就得到了均方根误差（Root Mean Squared Error）：

$$ RMSE：sqrt{frac{1}{m}sum_{i=1}^{m}(y_{test}^{(i)} - hat{y}\_{test}^{(i)})^2} = sqrt{MSE_{test}} $$

import math

math.sqrt(np.sum((y_predict - y_true) ** 2) / len(y_true))

还有一种衡量方法是平均绝对误差（Mean Absolute Error），对测试数据集中预测值与真值的差的绝对值取和，再取一个平均值：

$$ MAE：frac{1}{m}sum_{i=1}^{m}|y_{test}^{(i)} - hat{y}_{test}^{(i)}| $$

np.sum(np.absolute(y_predict - y_true)) / len(y_true)

注：Scikit Learn 的 metrics 模块中的 mean_squared_error() 方法表示 MSE，mean_absolute_error() 方法表示 MAE，没有表示 RMSE 的方法。

更近一步，MSE、RMSE 和 MAE 的局限性在于对模型的衡量只能做到数值越小表示模型越好，而通常对模型的衡量使用1表示最好，0表示最差，因此引入了新的指标：R Squared，计算公式为：

$$ R^2 = 1 - frac{SS_{residual}}{SS_{total}} $$

$SS_{residual} = sum_{i=1}^{m}(hat{y}^{(i)} - y^{(i)})^2$，表示使用模型产生的错误；$SS_{total} = sum_{i=1}^{m}(overline{y} - y^{(i)})^2$，表示使用 $y = overline{y}$ 预测产生的错误。

更深入的讲，对于每一个预测样本的 x 的预测值都为样本的均值 $overline{y}$ ，这样的模型称为基准模型；当我们的模型等于基准模型时，$R^2$ 的值为0，当我们的模型不犯任何错误时 $R^2$ 得到最大值1。

$R^2$ 还可以进行转换，转换结果为：

$$ R^2 = 1 - frac{MSE(hat{y}, y)}{Var(y)} $$

实现也很简单：

1 - np.sum((y_predict - y_true) ** 2) / len(y_true) / np.var(y_true)

注：Scikit Learn 的 metrics 模块中的 r2_score() 方法表示 R Squared。

当有不只一个自变量时，即为多元线性回归，如图：

对于有 n 个自变量来说，我们想获得的线性模型为：

$$ y = heta_{0} + heta_{1}x_{1} + heta_{2}x_{2} + ... + heta_{n}x_{n} $$

根据简单线性回归的思路，我们的目标即为：

找到 $ heta_{0}$，$ heta_{1}$，$ heta_{2}$，...，$ heta_{n}$，使得 $sum_{i=1}^{m}(y^{(i)} - hat{y}^2)^2$ 尽可能的小，其中 $hat{y}^{(i)} = heta_{0} + heta_{1}X_{1}^{(i)} + heta_{2}X_{2}^{(i)} + ... + heta_{n}X_{n}^{(i)}$ 。

$hat{y}^{(i)}$：训练数据中第 i 个样本的预测值；$X_{j}^{(i)}$：训练数据中第 i 个样本的第 j 个自变量。

如果用矩阵表示即为：

$$ hat{y}^{(i)} = X^{(i)}· heta $$

其中：${X^{(i)} = (X_{0}^{(i)},X_{1}^{(i)},X_{2}^{(i)},...,X_{n}^{(i)}), X_{0}^{(i)}equiv1}$；$ heta = ( heta_{0}, heta_{1}, heta_{2},..., heta_{n})^T$。

更进一步，将 $hat{y}^{(i)}$ 也使用矩阵表示，即为：

$$ hat{y} = X_b· heta $$

其中：$X_b = egin{pmatrix} 1 & X_1^{(1)} & X_2^{(1)} & cdots & X_n^{(1)} 1 & X_1^{(2)} & X_2^{(2)} & cdots & X_n^{(2)} cdots & & & & cdots 1 & X_1^{(m)} & X_2^{(m)} & cdots & X_n^{(m)} end{pmatrix}$，$ heta = egin{pmatrix} heta_0 heta_1 heta_2 cdots heta_n end{pmatrix} $

因此，我们目标就成了：使 $(y-X_b· heta)^T(y-X_b· heta)$ 尽可能小。而对于这个公式的解，称为多元线性回归的正规方程解（Nomal Equation）：

$$ heta = (X_b^TXb)^{-1}(X_b^Ty) $$

将多元线性回归实现在 LinearRegression 类中，且使用 Scikit Learn 的风格。

_init_() 方法首先初始化线性回归模型，_theta 表示 $ heta$，interception_ 表示截距，chef_ 表示回归模型中自变量的系数：

class LinearRegression:
    def __init__(self):
        self.coef_ = None
        self.interceiption_ = None
        self._theta = None

fit_normal() 方法根据训练数据集训练模型，X_b 表示添加了 $X_{0}^{(i)}equiv1$ 的样本特征数据，并且使用多元线性回归的正规方程解求出 $ heta$：

def fit_normal(self, X_train, 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

predict() 方法为预测方法，同样使用了矩阵乘：

def predict(self, X_predict):
    X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
    return X_b.dot(self._theta)

score() 根据给定的测试数据集使用 R Squared 指标计算模型的准确度：

def score(self, X_test, y_test):
    y_predict = self.predict(X_test)
    return r2_score(y_test, y_predict)

Scikit Learn 中的线性回归实现放在 linear_model 模块中，使用方法如下：
from sklearn.linear_model import LinearRegression

线性回归算法是典型的参数学习的算法，只能解决回归问题，其对数据具有强解释性。

缺点是多元线性回归的正规方程解 $ heta = (X_b^TXb)^{-1}(X_b^Ty)$ 的时间复杂度高，为 $O(n^3)$，可优化为 $O(n^{2.4})$。

Github | ML-Algorithms-Action