线性回归-最小二乘法

1、简单线性回归

所要求得的结果是一个具体的数值，而不是一个类别的话，则该问题是回归问题。只有一个特征的回归问题，成为简单线性回归。两个变量之间存在一次方函数关系，就称它们之间存在线性关系。拟合就是把平面上一系列的点，用一条光滑的曲线连接起来。
使用方程$y=ax+b$ 最大程度拟合样本点。图中每一个点对应了一个样本特征$x^i$，表示第$i$个数据点,样本对应的结果值为$y^i$。
如找到了拟合程度最大的方程，就是直到a，b的值。将$x^i$带入找到的这个方程，对应的结果值为：${\hat{y}}^i$。即：${\hat{y}}^{(i)}=a{x}^{(i)}+b$。
找到这个最佳的a和b，等价于使$\left|y^{(i)}-{\hat{y}}^{(i)}\right|$最小，即，使${\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2$最小，等价于使${\sum }_{i=1}^m{\left(y^i-{\hat{y}}^i\right)}^2$最小，等价于使${\sum }_{i=1}^m{\left(y^i-a{x}^i-b\right)}^2$最小。
这里我们称${\sum }_{i=1}^m{\left(y^i-a{x}^i-b\right)}^2$这个函数为损失函数。这里的损失可以理解为误差或错误。损失越小也就是错误越小。

2、使用最小二乘法求a，b

令损失函数为：$J(a,b)={\sum }_{i=1}^m{\left(y^{(i)}-a{x}^{(i)}-b\right)}^2$
问题转化为了求损失函数的最值问题。高等数学的求偏导：

令：$\frac{\partial J(a,b)}{\partial a}=0$ 和$\frac{\partial J(a,b)}{\partial b}=0$

对b求偏导：
因为：$\frac{\partial J(a,b)}{\partial b}$=$2{\sum }_{i=1}^m\left(y^i-a{x}^i-b\right)(-1)$

所以：$2{\sum }_{i=1}^m\left(y^i-a{x}^i-b\right)(-1)=0$

即：${\sum }_{i=1}^m\left(y^i-a{x}^i-b\right)=0$

即：${\sum }_{i=1}^m{y}^i-{\sum }_{i=1}^{m}a{x}^i-mb=0$

即：$mb={\sum }_{i=1}^m{y}^i-{\sum }_{i=1}^{m}a{x}^i$

即：$b=\bar{y}-a\bar{x}$
对a求偏导：
因为：$\frac{\partial J(a,b)}{\partial a}$=$2{\sum }_{i=1}^m\left(y^i-a{x}^i-b\right)\left(-x^i\right)$

所以：${\sum }_{i=1}^m\left(y^{(i)}-a{x}^{(i)}-b\right)x^{(i)}=0$

即：${\sum }_{i=1}^m\left(y^{(i)}-a{x}^{(i)}-\bar{y}+a\bar{x}\right)x^{(i)}=0$

此时上式仅有一个未知数a。整理后有：

${\sum }_{i=1}^m\left(x^{(i)}{y}^{(i)}-x^{(i)}\bar{y}-a{\left(x^{(i)}\right)}^2+a\bar{x}{x}^{(i)}\right)=0$

${\sum }_{i=1}^m\left(x^{(i)}{y}^{(i)}-x^{(i)}\bar{y}\right)-{\sum }_{i=1}^m\left(a{\left(x^{(i)}\right)}^2-a\bar{x}{x}^{(i)}\right)=0$

即：
$a=\frac{{\sum }_{i=1}^m\left(x^i{y}^i-x^i\bar{y}\right)}{{\sum }_{i=1}^m\left({\left(x^i\right)}^2-\bar{x}{x}^i\right)}$

其中：${\sum }_{i=1}^m{x}^i\bar{y}=\bar{y}{\sum }_{i=1}^m{x}^i=m\bar{y}\frac{1}{m}{\sum }_{i=1}^m{x}^i=m\bar{y}\cdot \bar{x}$

则：${\sum }_{i=1}^m{x}^i\bar{y}=m\bar{y}\cdot \bar{x}=\bar{x}{\sum }_{i=1}^m{y}^i={\sum }_{i=1}^m\bar{x}{y}^i$

即：${\sum }_{i=1}^m{x}^i\bar{y}={\sum }_{i=1}^m\bar{x}{y}^i$

那么：

$\begin{array}{rl}a& =\frac{{\sum }_{i=1}^m\left(x^i{y}^i-x^i\bar{y}\right)}{{\sum }_{i=1}^m\left({\left(x^i\right)}^2-\bar{x}{x}^i\right)}\\ & =\frac{{\sum }_{i=1}^m\left(x^i{y}^i-x^i\bar{y}-\bar{x}{y}^i+\bar{x}\cdot \bar{y}\right)}{{\sum }_{i=1}^m\left({\left(x^i\right)}^2-\bar{x}{x}^i-\bar{x}{x}^i+{\bar{x}}^2\right)}\end{array}$
则：
$\begin{array}{rl}a& =\frac{{\sum }_{i=1}^m\left(x^i{y}^i-x^2\bar{y}\right)}{{\sum }_{i=1}^m\left({\left(x^i\right)}^2-\bar{x}{x}^i\right)}\\ & =\frac{{\sum }_{i=1}^m\left(x^i{y}^i-x^i\bar{y}-\bar{x}{y}^i+\bar{x}\cdot \bar{y}\right)}{{\sum }_{i=1}^m\left({\left(x^i\right)}^2-2\bar{x}{x}^i+{\bar{x}}^2\right)}\\ & =\frac{{\sum }_{i=1}^m\left(x^i-\bar{x}\right)\left(y^i-\bar{y}\right)}{{\left(x^i-\bar{x}\right)}^2}\end{array}$

那么问题更加明确了：找到a、b的取值，使损失函数取最小值。这里：

$a=\frac{{\sum }_{i=1}^m\left(x^i-\bar{x}\right)\left(y^i-\bar{y}\right)}{{\sum }_{i=1}^m{\left(x^i-\bar{x}\right)}^2}$

$b=\bar{y}-a\bar{x}$
代码实现：

# 最小二乘法实现的
import numpy as np

class SimpleLinearRegression1:
   def __init__(self):
       self.a_ = None
       self.b_ = None
   
   def fit(self,x_train,y_train):
       x_mean = np.mean(x_train)
       y_mean = np.mean(y_train)
       
       num = 0.0
       d = 0.0
       # 先通过循环来求解
       for x,y in zip (x_train,y_train):
           num += (x - x_mean) * (y - y_mean)
           d += (x - x_mean) ** 2
       self.a_ = num / d
       self.b_ = y_mean - self.a_ * x_mean
       
       return self
   
   def predict(self,x_predict):
       return np.array([self._predict(x) for x in x_predict])
   
   def _predict(self,x_single):
       return self.a_ * x_single + self.b_
   
   def __repr__(self):
       return "SimpleLinearRegression1(a=%s,b=%s)" % (self.a_,self.b_)

其中加入for循环来计算效率太低了，其实可以使用numpy提供的向量运算来实现。

2.1 向量化运算

代码实现：

# 最小二乘法实现的
import numpy as np

class SimpleLinearRegression2:
    def __init__(self):
        self.a_ = None
        self.b_ = None
    
    def fit(self,x_train,y_train):
        x_mean = np.mean(x_train)
        y_mean = np.mean(y_train)
        
        # 通过向量化运算
        num = (x_train - x_mean).dot(y_train - y_mean) # 返回的是标量
        d = (x_train - x_mean).dot(x_train - x_mean)
        self.a_ = num / d
        self.b_ = y_mean - self.a_ * x_mean
        
        return self
    
    def predict(self,x_predict):
        return np.array([self._predict(x) for x in x_predict])
    
    def _predict(self,x_single):
        return self.a_ * x_single + self.b_
    
    def __repr__(self):
        return "SimpleLinearRegression2(a=%s,b=%s)" % (self.a_,self.b_)

3、评估线性回归的指标

在分类问题中可以使用准确度衡量，那在回归问题中使用什么衡量呢。
将数据分为训练数据和测试数据。利用训练集得到的a，b，将测试集的样本数据带入训练模型，也就是带入拟合后的方程，得到一个结果，对比模型得出的结果与真是的结果的差值。
在训练的时候：使${\sum }_{i=1}^m{\left(y_{\text{train}}^{(i)}-a{x}_{\text{train}}^{(i)}-b\right)}^2$尽可能小。训练结束后，得到a，b。即：得到相应的${\hat{y}}_{\text{test}}^{(i)}=a{x}_{\text{test}}^{(i)}+b$

3.1、均方误差(MSE,Mean Squared Error)

$\frac{1}{m}{\sum }_{i=1}^m{\left(y_{\text{test }}^{(i)}-{\hat{y}}_{\text{test }}^{(i)}\right)}^2$

3.2、均方根误差(RMSE,Root Mean Squared Error)

$\sqrt{\frac{1}{m}{\sum }_{i=1}^m{\left(y_{\text{test }}^{(i)}-{\hat{y}}_{\text{test }}^{(i)}\right)}^2}=\sqrt{MS{E}_{\text{test }}}$

3.3、平均绝对误差(MAE,Mean Absolute Error)

$\frac{1}{m}{\sum }_{i=1}^m\mid y_{\text{test}}^{(i)}-{\hat{y}}_{\text{test}}^{(i)}$

3.4 判定系数 R-square

上面我们介绍了MSE,RMSE和MAE这三个指标。其实这三个指标还有一定的问题。我们评价分类问题的指标是分类的准确度：1表示最好，0表示最差。即使我们的分类问题不同，我们也可以用这个指标来表示不同分类算分的优劣。但是上面的这三个指标是没有这个性质的。这个问题的解决方法是使用判定系数R Square)这个指标。
对于分子来说，我们可以理解为使用我们的模型预测产生的错误；对于分母来说，可以理解为使用均值预测产生的错误。如果一个模型使用均值来进行预测，那么这个模型叫基准模型(Baseline Model)，就是说我们的模型至少要比基准模型要好。如果整个式子的结果大于1，说明我们的模型甚至比不上基准模型。 结合整个式子来考虑的话， $R^2$越大越好，最大为 1 ，因为最好的模型就是使得分子为 0 ，整个式子为 1 。
这个式子还可以这样变形：

$R^2=1-\frac{{\sum }_i{\left({\hat{y}}^i-y^i\right)}^2}{{\sum }_i{\left(\bar{y}-y^i\right)}^2}=1-\frac{\left({\sum }_{i=1}^m{\left({\hat{y}}^i-y^i\right)}^2\right)/m}{\left({\sum }_{i=1}^m{\left(y^i-\bar{y}\right)}^2\right)/m}$
则： $R^2=1-MSE(\hat{y},y)/Var(y)$,分子变成了MSE，分母变成了y的方差。
代码实现：

def mean_squared_error(y_true,y_predict):
    return  np.sum((y_true - y_predict)**2)/ len(y_true)

def root_mean_squared_error(y_true,y_predict):
    return  sqrt(mean_squared_error(y_true,y_predict))

def mean_absolute_error(y_true,y_predict):
    return np.sum(np.absolute(y_true - y_predict)) / len(y_true)
    
def r_squire(y_test,y_predict):
	return 1-mean_squared_error(y_test,y_predict) / np.var(y_test)

最后，在自定义的简单线性回归中添加评分功能：

# 最小二乘法实现的
import numpy as np

class SimpleLinearRegression2:
    def __init__(self):
        self.a_ = None
        self.b_ = None
    
    def fit(self,x_train,y_train):
        x_mean = np.mean(x_train)
        y_mean = np.mean(y_train)
        
        # 通过向量化运算
        num = (x_train - x_mean).dot(y_train - y_mean) #dot是向量点乘，返回的是标量
        d = (x_train - x_mean).dot(x_train - x_mean)
        self.a_ = num / d
        self.b_ = y_mean - self.a_ * x_mean
        
        return self
    
    def predict(self,x_predict):
        return np.array([self._predict(x) for x in x_predict])
    
    def _predict(self,x_single):
        return self.a_ * x_single + self.b_
    
    def score(self,x_test,y_test):
        y_predict = self.predict(x_test)
        return mean_squared_error(y_test,y_predict) / np.var(y_test)
    
    def __repr__(self):
        return "SimpleLinearRegression2(a=%s,b=%s)" % (self.a_,self.b_)

4、多元线性回归

其中${\theta }_0$是截距，其余${\theta }_i$为各个系数。

问题转化为：求得一个$\theta$，使目标尽可能小。下边的式子称为多元线性回归的正规方程解(Normal Equation)

代码实现：

import numpy as np
from sklearn.metrics import r2_score


class LinearRegression:
    def __init(self):
        self.coef_ = None  # 系数
        self.interception_ = None  # 截距
        self._theta = None

    def fit_normal(self, X_train, y_train):
        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])  # 构造X_b X_train加上 虚构的都等于1的列
        self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)  # 通过正规方程解求得theta

        # 分开保存
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]
        return self


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

    def score(self, X_test, y_test):
        y_predict = self.predict(X_test)
        return r2_score(y_test, y_predict)
    
    def __repr__(self):
        return "LinearRegression(coef_=%s,interception_=%s)" % (self.coef_, self.interception_)

5、总结

线性回归算法只能解决回归问题，它是一个典型的参数学习问题。线性回归算法的特点就是对数据有很强的解释性。
我们在使用线性回归算法时，是有一个假设的。假设就是数据和最终的输出结果之间有一定的线性关系。
通过正规方程解来实现的线性回归算法， 如果我们的样本数量或特征数量巨大的话，使用这种方法是不合适的。梯度下降法是求解最优模型的通用方法。

来源：机器学习入门
最小二乘法的本质是什么