机器学习入门+实战初级（一）—— 线性回归

欢迎关注公众号K的笔记阅读博主更多优质学习内容

线性回归的概念及数学推导

一般假定

首先我们假定机器学习数据集的一般形式为： { ( x 1 , y 1 ) , ( x 2 , y 2 ) , … , ( x n , y n ) } \left\{\left(x_{1}, y_{1}\right),\left(x_{2}, y_{2}\right), \ldots,\left(x_{n}, y_{n}\right)\right\} {(x1​,y1​),(x2​,y2​),…,(xn​,yn​)}，其中$x_i=\left(x_{i1};x_{i2};x_{i3};\dots ;x_{id}\right),y_i\in R$，可以用以下形式描述 $x_i$ 与 $y_i$ 之间的关系：
$$\begin{array}{rlr}f(x)& ={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\dots +{\theta }_d{x}_d& =\sum _{i=0}^d{\theta }_i{x}_i\end{array}$$那么我们应该如何确定各个 $\theta$ 的值，使得 $f(x)$ 尽可能地与 $y_i$ 接近呢？一般我们会选择使得均方误差最小化的 $\theta$ 组合，即$$argmi{n}_{\theta }J(\theta )=\frac{1}{2}\sum _{j=1}^n{\left(h_{\theta }\left(x^{(i)}\right)-y^{(i)}\right)}^2$$
那么问题来了，这个一般是从何而来，我们为何要选择均方误差最小化，绝对值最小化等其他距离表示形式它不香吗？

均方误差最小化的数学推导

下面我们使用极大似然估计来推导我们为何要使用均方误差最小化这个指标。我们首先可以将目标值与变量写成如下形式：
$$y^{(i)}={\theta }^T{x}^{(i)}+{ϵ}^{(i)}$$
其中 $ϵ$ 为随即噪音，我们一般假定它独立同分布，服从高斯分布（由中心极限定理），由此我们有：$$p\left({ϵ}^{(i)}\right)=\frac{1}{\sqrt{2\pi }\sigma }exp\left(-\frac{{\left({ϵ}^{(i)}\right)}^2}{2{\sigma }^2}\right)$$也即为：$$p\left(y^{(i)}|x^{(i)};\theta \right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(y^{(i)}-{\theta }^T{x}^{(i)}\right)}^2}{2{\sigma }^2}\right)$$由于我们假定样本独立同分布，因此对于全样本，我们有：$$L(\theta )=p(\vec{y}|X;\theta )=\prod _{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }exp\left(-\frac{{\left(y^{(i)}-{\theta }^T{x}^{(i)}\right)}^2}{2{\sigma }^2}\right)$$在这里我们选择一个合适的 $\theta$ 使得似然函数最大化，这即是我们极大似然估计在做的事。
为了方便计算，我们计算时通常对对数似然函数求最大值：
$\begin{array}{rl}l(\theta )& =\log L(\theta )=\log \prod _{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(y^{(i)}-{\theta }^T{x}^{(i)}\right)}^2}{2{\sigma }^2}\right)\\ & =\sum _{i=1}^n\log \frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(y^{(i)}-{\theta }^T{x}^{(i)}\right)}^2}{2{\sigma }^2}\right)\\ & =\mathrm{nlog}\frac{1}{\sqrt{2\pi \sigma }}-\frac{1}{{\sigma }^2}\cdot \frac{1}{2}\sum _{i=1}^n({\left(y^{(i)}-{\theta }^T{x}^{(i)}\right)}^2\end{array}$
显然，最大化 $l(\theta )$ 即最小化 $\frac{1}{2}{\sum }_{i=1}^n({\left(y^{(i)}-{\theta }^T{x}^{(i)}\right)}^2$

线性回归损失函数、代价函数、目标函数

损失函数(Loss Function)：度量单样本预测的错误程度，损失函数值越小，模型就越好。
代价函数(Cost Function)：度量全部样本集的平均误差。
目标函数(Object Function)：代价函数和正则化函数，最终要优化的函数。

我们记 $$J(\theta )=\underset{f\in F}{\min }\frac{1}{n}\sum _{i=1}^{n}L\left(y_i,f\left(x_i\right)\right)+\lambda J(F)$$其中 $F$ 为所有可能的 $f$ 的集合，我们可以通过一系列的优化方法对其进行优化，比如梯度下降法：$${\theta }_j:={\theta }_j-\alpha \frac{\partial J(\theta )}{\partial \theta }$$即 $${\theta }_j={\theta }_j+\alpha \sum _{i=1}^n\left(y^{(i)}-f_{\theta }(x{)}^{(i)}\right)x_j^{(i)}$$
此方法称为批梯度下降法，我们当然也有选取当个点进行迭代的下降法，称为随机梯度下降法，即：
$$\theta =\theta +\alpha \left(y^{(i)}-f_{\theta }(x{)}^{(i)}\right)x^{(i)}$$这种方法的好处是当数据点很多时，该方法运行效率更高；缺点是，因为每次只针对一个样本更新参数，未必找到最快路径达到最优值，甚至有时候会出现参数在最小值附近徘徊而不是立即收敛。但当数据量很大的时候，随机梯度下降法经常优于批梯度下降法。
梯度下降法也有一定的局限性，只有当目标函数为凸函数时才能找到全局最小值，否则有可能只找到局部最优值。
另外，当数据量较小的时候，我们可以直接考虑矩阵求解，也即令：
$$J(\theta )=\frac{1}{2}(X\theta -Y{)}^T(X\theta -Y)$$的导数为0，容易解得：$\theta ={\left(X^{T}X\right)}^{(-1)}{X}^{T}Y$，但这种方法往往只适用于数据量较小的情况，当数据量较大时矩阵求逆成为一件难事。
根据最优化理论，我们还有牛顿法、拟牛顿法等方法对参数 $\theta$ 进行求解，在此不一一阐述。

线性回归的评价指标

均方误差(MSE): $\frac{1}{m}{\sum }_{i=1}^m{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2$
均方根误差(RMSE)： $\sqrt{MSE}=\sqrt{\frac{1}{m}{\sum }_{i=1}^m{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2}$
平均绝对误差(MAE): $\frac{1}{m}{\sum }_{i=1}^m|y^{(i)}-{\hat{y}}^{(i)}|.$
事实上更常用的是 $R^2$，其定义为
$$R^2:=1-\frac{{\sum }_{i=1}^m{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2}{{\sum }_{i=1}^m{\left(\bar{y}-{\hat{y}}^{(i)}\right)}^2}=1-\frac{\frac{1}{m}{\sum }_{i=1}^m{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2}{\frac{1}{m}{\sum }_{i=1}^m{\left(\bar{y}-{\hat{y}}^{(i)}\right)}^2}=1-\frac{MSE}{VAR}$$
我们可以将之理解为，回归模型可以成功解释的数据方差部分在数据固有方差中所占的比例，$R^2$ 越接近1，表示可解释力度越大，模型拟合的效果越好。

线性回归的 python 实例

首先我们要生成一个近似线性的数据集，代码如下：

import numpy as np
np.random.seed(1)
x = np.random.rand(500, 3)
y = x.dot(np.array([4.2,5.7,10.8]))

下面这段代码直接使用 LinearRegression 库进行线性回归：

from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
%matplotlib inline

lr = LinearRegression(fit_intercept=True)
lr.fit(x, y)
print('R2:%s' %(lr.score(x,y)))
x_test = np.array([2,4,5]).reshape(1,-1)
y_hat = lr.predict(x_test)
print("预测值为: %s" %(y_hat))

除此之外也可以采用矩阵求解：

class LR_LS():
	def __init__(self):
		self.w = None      
    def fit(self, X, y):
        self.w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
        def predict(self, X):
        y_pred = X.dot(self.w)
        return y_pred

if __name__ == "__main__":
    lr_ls = LR_LS()
    lr_ls.fit(x,y)
    print("估计的参数值：%s" %(lr_ls.w))
    x_test = np.array([2,4,5]).reshape(1,-1)
    print("预测值为: %s" %(lr_ls.predict(x_test)))

以及梯度下降法求解：

class LR_GD():
    def __init__(self):
        self.w = None     
    def fit(self,X,y,alpha=0.02,loss = 1e-10):
        y = y.reshape(-1,1) 
        [m,d] = np.shape(X) 
        self.w = np.zeros((d)) 
        tol = 1e5
        while tol > loss:
            h_f = X.dot(self.w).reshape(-1,1) 
            theta = self.w + alpha*np.mean(X*(y - h_f),axis=0)
            tol = np.sum(np.abs(theta - self.w))
            self.w = theta
    def predict(self, X):
        y_pred = X.dot(self.w)
        return y_pred  

if __name__ == "__main__":
    lr_gd = LR_D()
    lr_gd.fit(x,y)
    print("估计的参数值为：%s" %(lr_gd.w))
    x_test = np.array([2,4,5]).reshape(1,-1)
    print("预测值为：%s" %(lr_gd.predict(x_test)))

三种方法最后得到的结果应基本相同