【Machine Learning】回归分析之线性回归算法OLS

「爱情、让人受尽委屈。」 2021-10-30 01:34 405阅读 0赞

### 文章目录 ###

*  前言
 *  线性回归之回归算法OLS
 *   *  一、假设/原理
     *  二、经典例子
     *  三、建立模型
     *  四、完整求解思路
     *   *  4.1 求解误差
         *  4.2 误差分布假定
         *  4.3 似然函数求权重参数
         *   *  4.3.1 似然函数
             *  4.3.2 似然估计本质
             *  4.3.3 极大似然估计
             *  4.3.4 求导法获取最小二乘法的极小值
             *  4.3.4 小结
             *  4.3.5 梯度下降求权重参数
     *  五、源码实现
     *   *  5.1 前言
         *  5.2 算法步骤
     *  六、参考链接

# 前言 #

上篇博文我们介绍了[回归分析的一些入门的知识点][Link 1]，此篇博文我们将从最基本很典型的线性回归算法之OLS开始了解。

# 线性回归之回归算法OLS #

## 一、假设/原理 ##

## 二、经典例子 ##

我们要预测房屋的价值，其中考虑的特征包括房屋的面积，和房屋的已使用年限。如下图所示的4条样本：

现在，一个房屋面积为80，使用年限为5年的房屋，根据上表提示的数据预测下这个房屋的价值，这是我们的目标。  
思考：房屋面积和使用年限都会影响房屋的价值，不过我们现在还不知道它们各自对价值有多大的影响？此时我们预测的房屋价值是一个连续值，因此回归得到的是一个值，这是一个典型的二元回归问题，如果要从线性回归入手，就是二元线性回归。通俗点说就是找到一个面 ( x 1 , x 2 ) (x1, x2) (x1,x2)能很好的拟合（ y y y房屋价值）以上4个样本。

## 三、建立模型 ##

先从最简单的线性回归思路出发，这也是机器学习的基本思路，从最简单的模型入手。假设 θ 1 \\theta\_\{1\} θ1是房屋面积的权重参数， θ 2 \\theta\_\{2\} θ2是使用年限的权重参数，那么拟合的平面便可以表示为：  
 h θ ( x ) = ∑ i = 0 2 θ i x i h\_\{\\theta\}(x)=\\sum\_\{i=0\}^\{2\} \\theta\_\{i\} x\_\{i\} hθ(x)=i=0∑2θixi要习惯用矩阵的表达，上面这个求和公式用矩阵表达为：  
 θ T x \\theta^\{T\} x θTx  
其中 θ \\theta\_\{\} θ表示为：  
 ( θ 0 θ 1 θ 2 ) \\left(\\begin\{array\}\{l\}\{\\theta\_\{0\}\} \\\\ \{\\theta\_\{1\}\} \\\\ \{\\theta\_\{2\}\}\\end\{array\}\\right) ⎝⎛θ0θ1θ2⎠⎞  
 x x x表示为：  
 ( 1 x 1 x 2 ) \\left(\\begin\{array\}\{l\}\{1\} \\\\ \{x\_\{1\}\} \\\\ \{x\_\{2\}\}\\end\{array\}\\right) ⎝⎛1x1x2⎠⎞  
这里涉及到线性代数的相关矩阵转换的知识点，比较简单就不多说啦。

## 四、完整求解思路 ##

### 4.1 求解误差 ###

在假设了以上的模型后，接下来最重要的是求解方程中的3个参数，其中第一个参数为偏置项。我们知道预测值和真实值之间一般是存在误差的，误差值用 E \\mathcal\{E\} E表示，公式如下所示：  
 y ( i ) = θ T x ( i ) + ε ( i ) y^\{(i)\}=\\theta^\{T\} x^\{(i)\}+\\varepsilon^\{(i)\} y(i)=θTx(i)\+ε(i)写到这里有没有觉得此公式特别像一元线性函数呀，，，，  
其中： y ( i ) y^\{(i)\} y(i)是第 i i i个样本的真实值，注意这种标记方法，是很重要的奥。那么问题来了，误差的分布情况可以是任意的吗，还是需要满足某种分布规律才行？

### 4.2 误差分布假定 ###

以上这个问题是非常重要的，如果误差分布没有满足某个规律，这个就很难做出预测了，因为它没有规律啊！所以我们假设任何一个样本的误差项满足[独立同分布][Link 2] 并且服从均值为0，方差为一定值的[高斯分布][Link 2]。公式如下所示：  
 f ( x ) = 1 2 π σ exp ⁡ ( − ( x − μ ) 2 2 σ 2 ) f(x)=\\frac\{1\}\{\\sqrt\{2 \\pi\} \\sigma\} \\exp \\left(-\\frac\{(x-\\mu)^\{2\}\}\{2 \\sigma^\{2\}\}\\right) f(x)=2πσ1exp(−2σ2(x−μ)2) f ( x ) = 1 2 π σ e ( − x 2 2 σ 2 ) f(x)=\\frac\{1\}\{\\sqrt\{2 \\pi\} \\sigma\} e^\{\\left(-\\frac\{x^\{2\}\}\{2 \\sigma^\{2\}\}\\right)\} f(x)=2πσ1e(−2σ2x2)在做出这个假定，分布服从高斯分布后，我们就可以将误差项直接带入一维高斯分布的公式中。然后将误差项： y − θ T x y-\\theta^\{T\} x y−θTx 带入上式，可得：  
 f = 1 2 π σ e ( − ( y − θ T x ) 2 2 σ 2 ) f=\\frac\{1\}\{\\sqrt\{2 \\pi\} \\sigma\} e^\{\\left(-\\frac\{\\left(y-\\theta^\{T\} x\\right)^\{2\}\}\{2 \\sigma^\{2\}\}\\right)\} f=2πσ1e(−2σ2(y−θTx)2)上式中的 x x x和 y y y，方差都是已知量， f f f为概率的取值，那么，由这个公式该如何求解参数 θ \\theta θ呢？

### 4.3 似然函数求权重参数 ###

[似然函数][Link 3]的确是求解类似问题的常用解决方法，包括以后的解决其他模型的参数(权重参数)，也有可能用到[似然函数][Link 3]。

#### 4.3.1 似然函数 ####

首先构建似然函数 L ( θ ∣ x ) L(\\theta|x) L(θ∣x) ，假设一共有 m m m个房屋相关样本，那么进一步得到似然函数（它是参数 θ \\theta θ为自变量的函数，这个一定要注意了，似然函数将概率转化为似然，这个还是似然的强大之处了）：  
 L ( θ ) = ∏ i = 1 m 1 2 π σ exp ⁡ ( ( y ( i ) − θ T x ( i ) ) 2 2 σ 2 ) L(\\theta)=\\prod\_\{i=1\}^\{m\} \\frac\{1\}\{\\sqrt\{2 \\pi\} \\sigma\} \\exp ^\{\\left(\\frac\{\\left(y^\{(i)\}-\\theta^\{T\} x^\{(i)\}\\right)^\{2\}\}\{2 \\sigma^\{2\}\}\\right)\} L(θ)=i=1∏m2πσ1exp(2σ2(y(i)−θTx(i))2)上式的意思是 m 个样本的误差分布的概率乘积，这就是概率似然函数。  
提到似然函数，那不得不提最大似然函数估计吧，为什么呢？因为既然似然函数是关于误差分布的发生概率的乘积，既然这些  
分布值都已经实实在在的出现了，为什么不求出这种 θ \\theta θ ，它能使得事件尽可能地逼近样本值，这就是最大似然估计。

#### 4.3.2 似然估计本质 ####

本质便是根据已有的大量样本（实际上就是利用已知的条件）来推断事件本身的一些属性参数的方法，最大估计更是最能反映这些出现的样本的，所以这个参数值也是最可靠和让人信任的，得到这个参数值后，等来了一个新样本  X ( i + 1 ) X(i+1) X(i\+1)后，我们可以预测它的标签值。

#### 4.3.3 极大似然估计 ####

为了让上式最大，因为是各项相乘不好求最大值，想到取对数：称为对数似然，这样就转换为求和。  
 log ⁡ ( L ( θ ) ) \\log (L(\\theta)) log(L(θ))转化后的结果为：  
  const  − coeff ⁡ × ∑ i = 1 m ( y ( i ) − θ T x ( i ) ) 2 \\text \{ const \}-\\operatorname\{coeff\} \\times \\sum\_\{i=1\}^\{m\}\\left(y^\{(i)\}-\\theta^\{T\} x^\{(i)\}\\right)^\{2\}  const −coeff×i=1∑m(y(i)−θTx(i))2要想求这个式子的极大似然值，即极大值，也就是要求解coeff后的那项的极小值吧，就是下面这项：  
 J ( θ ) = ∑ i = 1 m ( y ( i ) − θ T x ( i ) ) 2 J(\\theta)=\\sum\_\{i=1\}^\{m\}\\left(y^\{(i)\}-\\theta^\{T\} x^\{(i)\}\\right)^\{2\} J(θ)=i=1∑m(y(i)−θTx(i))2上个式子有个很容易记得名字，叫做最小二乘项，现在清楚地推导出了最小二乘项 。

#### 4.3.4 求导法获取最小二乘法的极小值 ####

为了求解上个式子的极小值，首先想到的是求偏导 = 0 =0 =0，然后得出极小值。 J ( θ ) J(\\theta) J(θ)写成矩阵的形式：  
 J ( θ ) = ( X θ − y ) T ( X θ − y ) J(\\theta)=(X \\theta-y)^\{T\}(X \\theta-y) J(θ)=(Xθ−y)T(Xθ−y)其中 X X X在本文第三节中房屋那个例子中：4行3列的矩阵， θ \\theta^\{\} θ是不是3行1列的矩阵，他俩相乘是不是等于4行1列的矩阵，然后y是4行1列的矩阵，是不是正好能做减法啊。对其求导得到如下式子，中间的过程大家自行推导吧。  
 ∇ θ J ( θ ) = 2 ( X T X θ − X T y ) \\nabla\_\{\\theta\} J(\\theta)=2\\left(X^\{T\} X \\theta-X^\{T\} y\\right) ∇θJ(θ)=2(XTXθ−XTy)本小结第2个公式，J(theta)前有个1/2吗（我们只是没有写出来），这样前面的2实际上可以正好消除，也就是说J(theta)那个最小二乘项，最好带上那个1/2吧，也就是说：  
 J ( θ ) = 1 2 ∑ i = 1 m ( y ( i ) − θ T x ( i ) ) 2 J(\\theta)=\\frac\{1\}\{2\} \\sum\_\{i=1\}^\{m\}\\left(y^\{(i)\}-\\theta^\{T\} x^\{(i)\}\\right)^\{2\} J(θ)=21i=1∑m(y(i)−θTx(i))2再回头求完导后的那个式子： ∇ θ J ( θ ) = 2 ( X T X θ − X T y ) \\nabla\_\{\\theta\} J(\\theta)=2\\left(X^\{T\} X \\theta-X^\{T\} y\\right) ∇θJ(θ)=2(XTXθ−XTy)，令偏导等于0，不就是得到一个极小值吗，但是别忘了，这里又有一个假定： X T X X^\{T\} X XTX。假定这个为非奇异矩阵吧，因为只有非奇异矩阵才可求矩阵的逆啊！如果上面这项近似为奇异矩阵，那么就会引起一个最小二乘法的bug，这也是最小二乘法不能处理多重强相关性数据集的原因所在。  
假定不是奇异矩阵，那么参数theta这次可以求解出来了，即：  
 θ = ( X T X ) − 1 X T y \\theta=\\left(X^\{T\} X\\right)^\{-1\} X^\{T\} y θ=(XTX)−1XTy  
说明：具体的[最小二乘法矩阵求导过程][Link 4]

#### 4.3.4 小结 ####

在以上求解过程中做了一个 X T X X^\{T\} X XTX不能为奇异矩阵的假定，再加上之前的误差分布必须满足某种分布这个假定，所以最小二乘法直接求解得满足两个假定。以上我们通过数学的方法，借助似然函数，然后求似然函数对数的极大似然估计，直接把参数求出来了，这是必然？还是巧合？  
机器学习的参数一般是不能通过直接求解得出的，所以很明显是个巧合啊！那么如果不想用这种巧合的方法去求解，有没有更加通用的方法，来求解最小二乘项的极小值呢？答案是有的：梯度下降法。

#### 4.3.5 梯度下降求权重参数 ####

可参考博主的上一篇博文：[梯度下降算法介绍][Link 5]

## 五、源码实现 ##

### 5.1 前言 ###

最小二乘法适用于对处理的一堆数据，不必精确的经过每一点，而是根据图像到每个数据点的距离和最小确定函数。最小二乘法逼近的最简单的例子是根据一组观测值对 ( x 1 , y 1 ) , ( x 2 , y 2 ) … ( x n , y n ) (x1,y1),(x2,y2)…(xn,yn) (x1,y1),(x2,y2)…(xn,yn)来拟合一条直线。  
接下来我们就着重介绍一下一元线性回归的实现过程：  
1：公式  
 y = a x + b + ε y=a x+b+\\varepsilon y=ax\+b\+ε  
其中：  
 y y y为应变量(dependent variable), x x x自变量(independent variable) ， a a a为斜率( coeffient)， b b b为截距( intercept)， ε ε ε为误差，正态分布。线性回归的目的：找到一组 a 和 b a和b a和b，使得 ε ε ε最小。  
 y ^ = a x + b ε = y − y ^ \\begin\{array\}\{l\}\{\\hat\{y\}=a x+b\} \\\\ \{\\varepsilon=y-\\hat\{y\}\}\\end\{array\} y^=ax\+bε=y−y^ y ^ \\hat\{y\} y^读作y hat，也有人读作 y y y帽子。这里的帽子一般表示估计值，用来区别真实值 y y y。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0ppYWppa2FuZ19qams_size_16_color_FFFFFF_t_70]  
其中：  
金色的点是观测样本： y = a x + b + ε y=a x+b+\\varepsilon y=ax\+b\+ε  
红色的线为回归线： y ^ = a x + b \\hat\{y\}=a x+b\_\{\} y^=ax\+b  
绿色的线段为误差： ε = y − y ^ \\varepsilon=y-\\hat\{y\} ε=y−y^

### 5.2 算法步骤 ###

1.  a和b的起始值设置为零
2.  通过模型  y ^ = a x + b \\hat\{y\}=a x+b y^=ax\+b,就可计算出 y ^ \\hat\{y\} y^
3.  有了 y ^ \\hat\{y\} y^，就可以用优化方法算去更新参数
4.  重复2和3，直到找到J的最小值

""" author:jjk datetime:2019/5/2 coding:utf-8 project name:Pycharm_workstation Program function: https://mp.csdn.net/mdeditor/99714602 """
    
    ##最小二乘法
    import numpy as np  ##科学计算库
    import scipy as sp  ##在numpy基础上实现的部分算法库
    import matplotlib.pyplot as plt ##绘图库
    from scipy.optimize import leastsq ##引入最小二乘法算法
    
    ''' 设置样本数据，真实数据需要在这里处理 样本数据(Xi,Yi)，需要转换成数组(列表)形式 '''
    Xi=np.array([6.19,2.51,7.29,7.01,5.7,2.66,3.98,2.5,9.1,4.2])
    Yi=np.array([5.25,2.83,6.41,6.71,5.1,4.23,5.05,1.98,10.5,6.3])
    
    
    ''' 设定拟合函数和偏差函数 函数的形状确定过程： 1.先画样本图像 2.根据样本图像大致形状确定函数形式(直线、抛物线、正弦余弦等) '''
    
    ##需要拟合的函数func :指定函数的形状 说白了也就是模型
    def func(p,x):
      k,b=p
      return k*x+b
    # 或者
    def model(a,b,x):
        return a*x+b
    
    ##偏差函数：x,y都是列表:这里的x,y更上面的Xi,Yi中是一一对应的
    def error(p,x,y):
      return func(p,x)-y
    # 或者
    def cost_function(a,b,x,y):
        n = 5
        return 0.5/n * (np.square(y-a*x-b)).sum()
    
    
    ''' 主要部分：附带部分说明 1.leastsq函数的返回值tuple，第一个元素是求解结果，第二个是求解的代价值(个人理解) 2.官网的原话（第二个值）：Value of the cost function at the solution 3.实例：Para=>(array([ 0.61349535, 1.79409255]), 3) 4.返回值元组中第一个值的数量跟需要求解的参数的数量一致 '''
    
    #k,b的初始值，可以任意设定,经过几次试验，发现p0的值会影响cost的值：Para[1]
    p0=[1,20]
    
    #把error函数中除了p0以外的参数打包到args中(使用要求)
    Para=leastsq(error,p0,args=(Xi,Yi))
    
    #读取结果
    k,b=Para[0]
    print("k=",k,"b=",b)
    print("cost："+str(Para[1]))
    print("求解的拟合直线为:")
    print("y="+str(round(k,2))+"x+"+str(round(b,2)))
    
    ''' 绘图，看拟合效果. matplotlib默认不支持中文，label设置中文的话需要另行设置 如果报错，改成英文就可以 '''
    
    #画样本点
    plt.figure(figsize=(8,6)) ##指定图像比例： 8：6
    plt.scatter(Xi,Yi,color="green",label="样本数据",linewidth=2) #设置
    
    #画拟合直线
    x=np.linspace(0,12,100) ##在0-15直接画100个连续点
    y=k*x+b #函数式
    plt.plot(x,y,color="red",label="拟合直线",linewidth=2)
    plt.legend(loc='lower right') #绘制图例
    plt.show()

![在这里插入图片描述][20190819202517734.png]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0ppYWppa2FuZ19qams_size_16_color_FFFFFF_t_70 1]  
也可参考：[线性回归演示][Link 6]

## 六、参考链接 ##

1、[python机器学习手写算法系列——线性回归][python]  
2、[Python中实现最小二乘法思路及实现代码][Python]  
3、[Python实现基于最小二乘法的线性回归][Python 1]  
4、[机器学习-回归分析之线性回归][-]  
5、[机器学习算法-------线性回归法][-------]  
6、[最小二乘法-python实现][-python]  
7、[线性回归案例][Link 7]  
源码获取：  
1、[https://github.com/jiajikang1993/NLP\_related\_algorithm\_learned/tree/master/Machine%20Learning(ML)][https_github.com_jiajikang1993_NLP_related_algorithm_learned_tree_master_Machine_20Learning_ML]  
2、[一个下载数据集的网址—不一定是您需要的奥][Link 8]

[Link 1]: https://blog.csdn.net/Jiajikang_jjk/article/details/99720086
[Link 2]: https://blog.csdn.net/Jiajikang_jjk/article/details/99706496#_267
[Link 3]: https://blog.csdn.net/Jiajikang_jjk/article/details/99706496#_284
[Link 4]: https://blog.csdn.net/weixin_37606743/article/details/79111300
[Link 5]: https://blog.csdn.net/Jiajikang_jjk/article/details/99731119
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0ppYWppa2FuZ19qams_size_16_color_FFFFFF_t_70]: /images/20211029/2000b1d84b3746688a28580193052b39.png
[20190819202517734.png]: /images/20211029/6584aee99f8646c0b1a02681a874c5b8.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0ppYWppa2FuZ19qams_size_16_color_FFFFFF_t_70 1]: /images/20211029/280818f9f8cb4aabb05d4fc05dd52845.png
[Link 6]: https://github.com/juwikuang/machine_learning_step_by_step/blob/master/linear_regression.ipynb
[python]: https://blog.csdn.net/juwikuang/article/details/78420337
[Python]: https://www.jb51.net/article/132137.htm
[Python 1]: https://www.cnblogs.com/geeksongs/archive/2019/03/26/10602630.html
[-]: https://blog.csdn.net/Daycym/article/details/79922176
[-------]: https://blog.csdn.net/qq_34374664/article/details/80304588
[-python]: https://blog.csdn.net/deramer1/article/details/79055281
[Link 7]: https://blog.csdn.net/Daycym/article/details/80072528
[https_github.com_jiajikang1993_NLP_related_algorithm_learned_tree_master_Machine_20Learning_ML]: https://github.com/jiajikang1993/NLP_related_algorithm_learned/tree/master/Machine%20Learning%28ML%29
[Link 8]: http://archive.ics.uci.edu/ml/machine-learning-databases/