Linear Regression模型简介

一、概述

1、背景

• 最近几十年，人工智能经历了一轮又一轮的高潮和低谷；今天，机器学习、深度学习再一次被赋予强人工智能的历史使命；

• 机器学习作为人工智能中最重要的一环，天生就带有数学家和工程师的基因，了解机器学习的常见模型和算法还是很有意义的。

• 本文结合自己的学习情况介绍对线性回归模型的认识(如文章有不妥之处，望不吝赐教)

2、介绍

• 机器学习根据训练数据方法的不同，可以分为监督学习、半监督学习、无监督学习和强化学习。其中监督学习的主要任务是分类和回归

• 分类和回归问题都需要根据训练样本找到一个实数值函数g(x)；然后回归根据g(x)，可以推导出给定新样本的输出值，而分类是根据g(x)，推导出新样本的类别（如0或1）。

• 在现实世界中普遍存在变量之间的关系，关系大概分两种：确定的函数关系 或 统计相关关系（非确定但又相互依赖的关系）；回归分析就是研究变量之间统计相关关系的方法。

• 唯物辩证法中认为：事物之间的联系是普遍存在的、相互影响、相互作用和相互制约的；佛家认为：前世因，今世果；数学家认为： y = a_1x_1 + a_2x_2 + a_3x_3 + {\cdots}+ b

• 在回归问题中，按照自变量和变量之间的关系类型，分为线性回归分析和非线性回归分析（如对数模型:f(x) = \frac{1}{1+ e^{-x}})，本文介绍对象是线性回归（Linear Regression）模型。

3、线性回归模型

主要介绍：回归方程、损失函数、损失函数求解方法、过拟合和欠拟合、多项式

二、回归方程

1、概述

在回归分析中，只包括一个自变量和一个因变量，且两种的关系可以使用一条直线近似表示，称为一元线性回归分析。如果有两个及以上的自变量，且自变量和因变量是线性关系，称为多元线性回归分析。

• 一元线性回归方程： y = ax + b

• 多元线性回归方程： y = a_1x_1 + a_2x_2 + a_3x_3 + {\cdots} + b

2、建立回归方程

h_θ(x) = θ_0 + θ_1x_1+θ_2x_2+{\cdots}+θ_nx_n = \vec{θ}^T\vec{X}

其中，\vec{θ} = \begin{bmatrix} θ_0 \\ θ_1 \\ θ_2 \\ {\cdots} \\ θ_n \end{bmatrix} \,\,\,\,\,\,\vec{X} = \begin{bmatrix} 1 \\ x_1 \\ x_2 \\ {\cdots} \\ x_n \end{bmatrix}

说明：θ_0 , θ_1,θ_2,{\cdots}θ_n 是回归系数，我们的目标是求解出它们，然后代入回归方程，来预测新数据的输出值。

二、损失函数 （Loss Function）

1、公式

确立了回归方程后，我们使用均方误差（Mean Square Error,MSE）来评估真实值 与预测值之间的差异，公式如下：

J(θ) = \frac {1}{2n} \sum_{i=1}^n(h_θ(x^i) - y^i)^2,\,\,\text{n为样本数}

其矩阵表达式为：

J(θ) = \frac {1}{2n}(Xθ - Y)^T(Xθ - Y) ,\,\,\text{n为样本数}

说明：以上就是线性回归的损失函数。下面介绍公式的推导。

2、解释

• 假设特征的预测值和真实值的误差是\xi^{(i)} ,那么预测值θ^T x^{(i)}和真实值y^{(i)}满足如下公式：

y^{(i)} = θ^T x^{(i)} + \xi^{(i)}\,\,\,\,\,\,\,(1)

• 假设误差是独立的，随机的，无穷大的，根据中心极限定理，误差服从均值为\mu，方差为\sigma^2的正态分布，其概率分布如下:

p(\xi^{(i)}) = \frac {1}{\sqrt{2\pi\sigma}} exp(-\frac{{\xi^{(i)}}^2}{2\sigma^2})\,\,\,\,\,\,\,(2)

• 由公式(1)可知：\xi^{(i)} = y^{(i)} - θ^T x^{(i)}，代入公式(2)，得到一个y^{(i)}的概率公式如下：

p(y^{(i)} |x^{(i)};θ) = \frac {1}{\sqrt{2\pi} \sigma} exp(- \frac{(y^{(i)} - θ^T x^{(i)})^2}{2\sigma^2})\,\,\,\,\,\,\,(3)

• 对应的极大似然函数是：

L(θ) = \prod_{i=1}^np(y^{(i)} |x^{(i)};θ) = \prod_{i=1}^n\frac {1}{\sqrt{2\pi} \sigma} exp(- \frac{(y^{(i)} - θ^T x^{(i)})^2}{2\sigma^2})|x^{(i)};θ) \,\,\,\,\,\,\,(4)

• 极大似然函数取对数，求导后得到θ的极大似然估计如下

\hat{θ} = \frac {1}{2} \sum_{i=1}^n(h_θ(x^i) - y^i)^2 \,\,\,\,\,\,\,(5)

说明：极大似然函数估计知道，在\hat{θ} = \frac {1}{2} \sum_{i=1}^n(h_θ(x^i) - y^i)^2时，误差最小；为了后面计算推导，引入归一化的系数n, 最终得到的损失函数为：

J(θ) = \frac {1}{2n} \sum_{i=1}^n(h_θ(x^i) - y^i)^2,\,\,\text{n为样本数}

三、损失函数求解法1：普通最小二乘法(OLS)

1、求解

• 当矩阵X可逆时，可以用普通最小二乘法(OLS) ，求解出回归系数，公式如下(使用矩阵形式)：

\hat{θ} = (X^TX)^{-1}X^TY

• 对应的代码实现如下：

 def standRegre(xMat, yMat):
  '''
  :param xMat: X矩阵
  :param yMat: y矩阵
  :return: 回归系数矩阵
  '''

  xTx = xMat.T * xMat
  # linalg.det(xTx)计算行列式,行列式值为0.矩阵不可逆
  if linalg.det(xTx) == 0.0:
      print("This matrix is singular,cannot do inverse")
      return

  # .I用作求矩阵的逆矩阵
  ws = xTx.I * dot(xMat.T,yMat)
  return ws

• 回归方程对数据点的拟合情况如下：

效果图.png

说明：要求矩阵X必须满足可逆条件，否则是不可以解的。

2、解释

我们知道，损失函数的矩阵表达式为

J(θ) = \frac {1}{2n}(Xθ - Y)^T(Xθ - Y)
扩展后得到：
J(θ) = \frac {1}{2n}(θ^TX^TXθ - 2θ^TX^TY + YY )
对θ求偏导，令导数为0，得：
\frac{\partial J(θ)}{\partial θ} = \frac {1}{2n}(2θ^TX^TX - 2X^TY) = 0
最后求解得到：
θ^T = (X^TX)^{-1}(X^TY)

四、损失函数求解法2：梯度下降（Gradient Descent）算法

1、 简介

梯度下降算法是一种求局部最优解的方法；在数学上，函数的梯度(导数)方向是函数值增长最快的方向，其相反方向就是下降最快的方向。

• 首先对\vec{θ}赋值，这个值可以是随机的，也可以是一个全零的向量。

• 改变\vec{θ}的值，使得J(θ)按梯度下降的方向进行减少。

θ = θ - \alpha \frac{\partial J(θ)}{\partial θ}, \,\,\, \alpha\text{是学习率，步长}

   具体对θ的每个分量的求导公式如下：

\frac{\partial J(θ)}{\partial θ_j} = \frac{1}{n}\sum_{i=1}^n(h_θ(x^i) - y^i)x^i_j

说明：梯度下降（GD）是最小化损失函数的一种常用方法，主要有批量梯度下降和随机梯度下降两种迭代求解思路。

2、 批量梯度下降 (Batch Gradient Descent)

repeat until convergence {
\,\,\,\,\,\,\,\,\,θ_j:= θ_j - \alpha \frac{1}{n}\sum_{i=1}^n(h_θ(x^i) - y^i)x^i_j
}

说明：每更新一个θ_j ，都要遍历一次样本集；如果样本的数量n很大，开销比较大；BGD可以得到一个全局最优解，使得损失函数的风险最小

代码实现：

def bgd(rate, maxLoop, epsilon, X, Y):
    '''
    批量梯度下降法
    :param rate:  学习率
    :param maxLoop:
    :param epsilon: 收敛精度
    :param X: 样本特征矩阵
    :param Y: 标签矩阵
    :return: (theta, errors, thetas), timeConsumed
    '''

    # m是样本数，n是特征个数
    m, n = X.shape
    # 初始化theta,n*1维的零向量
    theta = np.zeros((n, 1))
    count = 0
    # 是否收敛
    converged = False
    error = float('inf')
    errors = []
    thetas = {}

    for j in range(n):
        thetas[j] = [theta[j, 0]]

    while count <= maxLoop:
        if (converged):
            break
        count = count + 1
        # sum = 0
        for j in range(n):
            # deriv = (X * theta - y).T * X[:, j] / m
            # print(deriv)
            deriv = np.float128(0)
            for i in range(m):
                deriv += ((X[i] * theta - Y[i]) * X[i,j])[0,0]/ m
            theta[j, 0] = theta[j, 0] - rate * deriv
            thetas[j].append(theta[j, 0])

        #误差
        error = loss_function(theta, X, y)
        errors.append(error[0, 0])
        # 如果已经收敛
        if (error < epsilon):
            converged = True
    return theta, errors, thetas

批量梯度下降.png

3、随机梯度下降 (Stochastic Gradient Descent)

repeat until convergence {
   for i=1 to n {
\,\,\,\,\,\,\,\,θ_j:= θ_j - \alpha (h_θ(x^i) - y^i)x^i_j
   }
}

说明1：在随机梯度下降法中，每更新θ_j 只需要一个样本 (x^i,y^i)；在样本数量巨大时，SGD的性能比较好，但是最终结果未必是全局最优解。

代码实现

def sgd(rate, maxLoop, epsilon, dataMat, labelMat):
  '''
  随机梯度下降法
  :param rate:  学习率
  :param maxLoop:
  :param epsilon: 收敛精度
  :param dataMat: 样本特征矩阵
  :param labelMat: 标签矩阵
  '''
  sampleCount, featureCount = dataMat.shape
  #初始化theta,featureCount * 1维的0向量
  theta = np.zeros((featureCount,1))
  count = 0
  converged = False
  error = float('inf')
  errors = []
  thetas = {}  #字典

  for j in range(featureCount):
      thetas[j] = [theta[j,0]]

  # 迭代到最大次数结束
  while count <= maxLoop:
      if (converged):
          break
      count = count + 1
      errors.append(float('inf'))
      # 获取样本
      for i in range(sampleCount):
          if (converged):
              break
          #计算实际和预测的误差
          diff = labelMat[i,0] - regression_function(theta,dataMat[i].T)
          # 通过第i个样本来优化theta
          for j in range(featureCount):
              theta[j, 0] = theta[j, 0] + rate * diff * dataMat[i,j]
              thetas[j].append(theta[j, 0])
          error = loss_function(theta, dataMat, labelMat)
          errors[-1] = error[0, 0]
          # 如果已经收敛
          if (error < epsilon):
              converged = True
   return theta, errors, thetas

随机梯度下降.png

4、其他

• 随机梯度下降法在还能够应用到一些更加复杂的算法中，如逻辑回归（Logic Regression）算法(虽有回归两字，但实际上是分类算法)等；

• 在回归模型的损失函数求解中，还有正规方程的方式，它不需要调节学习率\alpha，但是要求矩阵可逆，数据维度不能太大。

五、特征缩放

当样本的特征的数值之间差异比较大时，梯度下降的过程不仅曲折，而且耗时；这时候需要将各个特征归一化到统一的区间。常见的做法如下：

1、0均值标准化（Z-score normalization）

• 将原始数据集归一化为均值为0、方差为1的数据集，归一化后的特征将分布在 [−1,1] 区间，转化函数如下：

z=\frac {x - μ}{δ}

    其中，μ、σ分别为原始数据集的均值和标准差。

    说明：该种归一化方式要求原始数据的分布可以近似为高斯分布，否则归一化的效果会变得很糟糕。

• 代码实现如下：

def standardize(xArr):
    '''
    特征标准化处理
    :param xArr: 样本特征集
    :return: 标准化后的样本特征集
    '''
    m, n = xArr.shape
    # 归一化每一个特征
    for j in range(n):
        meanVal = xArr[:, j].mean(axis=0)
        std = xArr[:, j].std(axis=0)
        if std != 0:
            xArr[:, j] = (xArr[:, j] - meanVal) / pow(std,2)
        else:
            xArr[:, j] = 0
    return xArr

2、线性函数归一化（Min-Max Scaling）

• 实现对原始数据的等比例缩放，将原始数据线性化的方法转换到[0,1]的范围，公式如下：

X_{normal} = \frac {X - X_{min}}{ X_{max} - X_{min}}

    其中X_{norm}为归一化后的数据，X为原始数据，X_{max} 、X_{min}分别为原始数据集的最大值和最小值。

• 代码实现如下：

def normalize(xArr):
  """特征归一化处理
  Args:
      X: 样本特征集
  Returns:
      归一化后的样本特征集矩阵
  """
  m, n = xArr.shape
  # 归一化每一个特征
  for j in range(n):
      minVal = xArr[:, j].min(axis=0)
      maxVal = xArr[:, j].max(axis=0)
      diff = maxVal - minVal
      if diff != 0:
        xArr[:, j] = (xArr[:, j] - minVal) / diff
      else:
        xArr[:, j] = 0
   return xArr

3、总结

• 在分类、聚类算法中，需要使用距离来度量相似性的时候、或者使用PCA技术进行降维的时候，Z-score standardization表现更好。

• 在不涉及距离度量、协方差计算、数据不符合正太分布的时候，可以使用Min-Max Scaling或其他归一化方法。比如图像处理中，将RGB图像转换为灰度图像后将其值限定在[0 255]的范围。

六、多项式回归

1、多项式

多项式由称为不定元的变量和常数系数通过有限次加减法、乘法以及自然数幂次的乘方运算得到的代数表达式，可以分为一次一元多项式和多元多项式。我们主要讨论一元多项式，其形式如下：

y = a_0 + a_1 * x + a_2 * (x^2) + ... + a_n * (x^n) + e

2、多项式回归

• 当回归函数不能用直线来描述时，要考虑用非线性回归函数。 多项式回归属于非线性回归的一种。这里主要介绍一元多项式回归。

• 一般非线性回归函数是未知的，或即使已知也未必可以用一个简单的函数变换转化为线性模型。常用的做法是用因子的多项式。

• 从散点图观察到回归函数有一个“弯”，则可考虑用二次多项式；有两个弯则考虑用三次多项式；有三个弯则考虑用四次多项式，等等。

• 真实的回归函数未必就是某个次数的多项式，但只要拟合得好，用适当的多项式来近似真实的回归函数是可行的。

七、拟合的问题

1、概念

• 拟合：把平面上一系列的点，用一条光滑的曲线连接起来，拟合的曲线一般可以用函数表示，拟合曲线不是唯一的；在实践中，常常会出现过拟合和欠拟合问题；

欠拟合、拟合和过拟合.png

• 最左边的直线没有很好拟合数据，这类情况是欠拟合(underfitting)；最右边的曲线似乎对数据点做了很好的拟合(曲线通过了所有的点)，但是它无法泛化到新的数据样本中，是一个差的模型，这种情况是过拟合的，只有中间的曲线才是理想的拟合曲线。

• 泛化指的是一个假设模型能够应用到新样本的能力。新样本数据是指没有出现在训练集中的数据。

• 特征过多 或者 样本量少的情况都容易发生过拟合的情况。解决的方法有两种：1）减少特征个数，选择最重要的特征，抛弃一部分特征；2）正则化，保留所有的特征，但是会减小特征变量的数量级。