机器学习之线性回归

    在机器学习中，线性回归(Linear regression)是指利用最小二乘法对一个或多个自变量和因变量之间关系进行建模的一种回归分析方法。

线性回归初识

    在有监督的机器学习中，根据样本中$label$是离散值还是连续值，将学习问题分为分类问题和回归问题。线性回归主要是根据样本数据的线性加权拟合出一条直线，使得拟合的直线能逼近真实值。线性回归建模后得到的预测值用数学表达式表示为：

$$\bf{h_\theta(x)=\sum_{i=0}^n\theta_ix_i=\theta^T x}$$

注意：线性指的是参数是线性的，样本数据可以是非线性的。

线性回归目标函数

几个概念

中心极限定理

    实际问题中，很多随机现象可以看做是众多因素的独立影响的综合反映，往往近似服从正太分布。

最大似然估计

    最大似然估计就是利用已知的样本的结果，在使用某个模型的基础之上，反推导致这样结果的参数值。

最小二乘目标函数推导

    预测值和真实值之间关系表达式可以表示为如下（其中$\varepsilon^{(i)}$代表的是残差项）：

$$y^{(i)}=\theta^Tx^{(i)}+\varepsilon^{(i)}$$

假设数据集中有$m$个样本都是独立同分布的，即假设真实值只和自己样本数据有关，$\varepsilon^{(i)}$服从一定的分布，由于残差项产生的原因可能是由于建模过程中忽略了细微的特征，根据中心极限定理，残差项应该服从正常状态的分布：高斯分布（均值为0，方差为$\sigma^2$）
    已知$\varepsilon$服从正太分布，正太分布的参数是移植的，根据最大似然估计，反推$\theta$.
$\varepsilon$服从正太分布，用概率表示为：

$$p(\varepsilon^{(i)})=\frac{1}{\sqrt{2\pi \sigma}}exp(-\frac{(\varepsilon^{(i)})^2}{2\sigma^2})$$ $$p(y^{(i)}|x{(i)};\theta)=\frac{1}{\sqrt{2\pi \sigma}}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})$$

似然函数表示为：

$$L(\theta)=\prod_{i=1}^mp(y^{(i)}|x^{(i)};\theta)=\prod_{i=1}^m\frac{1}{\sqrt{2\pi \sigma}}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})$$

对数似然函数表示为：

$$l(\theta)=log L(\theta)=log\prod_{i=1}^mp(y^{(i)}|x^{(i)};\theta)$$ $$=log\prod_{i=1}^m\frac{1}{\sqrt{2\pi \sigma}}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})$$ $$=\sum_{i=0}^mlog\frac{1}{\sqrt{2\pi \sigma}}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})$$ $$=mlog\frac{1}{\sqrt{2\pi \sigma}}-\frac{1}{\sigma^2}·\frac12\sum_{i=1}^m(y^{(i)}-\theta^Tx^{(i)})^2$$

在上述数学式字当中，$\sigma$是常数，要想求得$l(\theta)$取得最大，则$\frac12\sum_{i=1}^m(y^{(i)}-\theta^Tx^{(i)})^2$需要取得最小值，因此得出最小二乘法的目标函数。
在线性回归中，可以得出损失函数，即所谓的目标函数是：

$$\bf{J(\theta)=\frac12\sum_{i=1}^m(y^{(i)}-\theta^Tx^{(i)})^2}$$

目标函数计算

线性回归目标函数：

$$\bf{J(\theta)=\frac12\sum_{i=1}^m(y^{(i)}-\theta^Tx^{(i)})^2}$$

用矩阵形式表示为：

$$\bf{J(\theta)=\frac12(X\theta-y)^T (X\theta-y)}$$

我们的目标是计算在$\theta$取何值时，目标函数$J(\theta)$取得最小值。对于凸函数来说，可以将求最小值问题转化为求驻点问题，对$\theta$求偏导，求等于目标函数的梯度：

$$\nabla_{\theta}J(\theta)=\nabla_{\theta}(\frac12(X\theta-y)^T(X\theta-y))$$ $$=\nabla_{\theta}(\frac12(\theta^TX^T-y^T)(X\theta-y))$$ $$=\nabla_{\theta}(\frac12(\theta^TX^TX\theta-\theta^TX^Ty-y^TX\theta+y^Ty))$$ $$=\frac12(2X^TX\theta-X^Ty-(y^TX)^T)=X^TX\theta-X^Ty$$

求驻点,若$X^TX$可逆，得：

$$\theta=(X^TX)^{(-1)}X^Ty$$

若$X^TX$不可逆或者是防止过拟合加入$\lambda$扰动可得：

$$\theta=(X^TX+\lambda I)^{(-1)}X^Ty$$

补充：
$X^TX$是半正定的，对于任意的非零向量$u$:

$$u^TX^TXu=(Xu)^TXu\ge0$$

加入扰动因子之后，$X^TX+\lambda I$是正定的，正定一定可逆。

修正目标函数:Ridge&LASSO

    有的时候，为了使目标函数尽可能的小，利用数据的高阶组合来达到这样的目的，但是这样导致的后果就是系数的值很大，过拟合，为了避免出现过拟合情况，会引入惩罚项因子$\lambda$。

Ridge

    Ridge是在优化过程的目标函数中使用L2 penalty,旨在把系数变小一些，但并非完全为0。通俗点理解，就是给原先的目标函数加一个约束条件，让各个系数的平方和小于某一个数（即规定范围为一个圆），求带有约束条件的极小值问题一般用拉格朗日乘子法。最终体现在数学表达式上，就是：

$$J(\theta)=\frac12\sum_{i=1}^m(h_{\theta}(x^{(i)})-y^{(i)})^2+\lambda\sum_{j=1}^n\theta_j^2$$

对目标函数求梯度，即可得到之前，加入扰动因子所得到的对$\theta$的估计：

$$\theta=(X^TX+\lambda I)^{(-1)}X^Ty$$

LASSO

    LASSO是在优化过程的目标函数中使用L1 penalty。LASSO不同于Ridge，它给原先求得的目标函数加的约束条件，是让各个系数的绝对值小于某一个数（范围是一个菱形），因为约束域有棱有角，所以有些系数会变为0，实现了特征选择的能力。其数学表达式可以为：

$$J(\theta)=\frac12\sum_{i=1}^m(h_{\theta}(x^{(i)})-y^{(i)})^2+\lambda \sum_{j=0}^m|\theta_j|$$

Elastic Net

    在实际工程上，Ridge在一些评价标准上优于LASSO，但是没有特征选择能力，因此，有人引入Elastic Net，用来折中LASSO和Ridge。数学表达式：

$$J(\theta)=\frac12\sum_{i=0}^m(h_{\theta}(x^{(i)})-y^{(i)})^2+\lambda(\rho\sum_{j=1}^m|\theta_j|+(1-\rho)\sum_{j=0}^m\theta_j^2)$$

问题1：采用L1 penalty，如何处理梯度呢?

参数和超参数

    参数($\theta$)估计：可以根据训练样本求得$\theta$，得到模型。
    超参数($\lambda$)选择：超参数就是参数的参数，给了超参数，参数可以求出来，超参数是无法根据样本求得的，人为可以选定候选集。根据超参数候选集，可以得出候选参数，进而可以在验证集上根据某种机器学习的评价指标，比如上一个博客所提到的SSE挑选效果较好的参数所对应的超参数，最终得到模型的参数和超参数。
    交叉验证：拿十折交叉验证来说，就是将训练集分为10份，九份用来训练，一份用来验证，每一份样本都会有一次当做验证数据，最终取平均值，来看效果。

梯度下降算法

    高维空间中的矩阵逆运算是复杂的，作为机器学习系列的博客，参数可以采用梯度下降算法来不断求得。

算法流程：

1. 随机初始化$\theta$
2. 沿着梯度方向(在某点增长最快的方向)的负方向不断迭代，使得更新后的目标函数$J(\theta)$更小，$\theta$的更新公式： $$\theta=\theta-\alpha·\frac{\partial J(\theta)}{\partial \theta}$$ 其中$\alpha$是学习率。
       采用梯度下降算法求1-100的平方根，点击代码下载。

   批量梯度下降

       在这个算法中，每次更新需要用到所有的训练样本。更新公式可以表示为： $$\theta_j=\theta_j-\alpha\sum_{i=0}^m(h_{\theta}(x^{(i)})-y^{(i)})x_j^{(i)}$$

   随机梯度下降

       在这个算法中，每次更新只用到一个训练样本。更新公式： $$\theta_j=\theta_j-\alpha(h_{\theta}(x^{(i)})-y^{(i)})x_j^{(i)}$$     随机梯度下降和批量梯度下降都有各自的优缺点，在训练速度方面，随机梯度下降更快，在准确度方面，随机梯度下降导致解很有可能不是最优。因此一个比较这种的方法是mini-batch梯度下降。

   mini-batch梯度下降

       在这个算法中，选用一部分样本的平均梯度作为更新方向。更新公式： $$\theta_j=\theta_j-\alpha\sum_{i=0}^t(h_{\theta}(x^{(i)})-y^{(i)})x_j^{(i)}$$

   LinearModel python版应用代码

   更多详细代码请参考个人github代码仓库
   https://github.com/MichelleYang2017/BlogCode

   LinearRegression 代码

   import pandas as pd
   import matplotlib as mpl
   import matplotlib.pyplot as plt
   from sklearn.model_selection import train_test_split
   from sklearn.linear_model import LinearRegression
   import  numpy as np
   
   if __name__=='__main__':
       filepath='./Data/Advertising.csv'
       data=pd.read_csv(filepath)#读取数据
       X=data[['TV','Radio','Newspaper']]
       y=data['Sales']
   
       #用来显示正常中文标签
       mpl.rcParams['font.sans-serif']=['simHei']
       #用来显示正常负号
       mpl.rcParams['axes.unicode_minus']=False
   
       
       X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.2,random_state=1)
       # print(X_train.shape,X_test.shape)
       '''
       fit_intercept是否有截距
       normalize:是否将数据归一化，减去均值/L2
       copy_X:当True时，X将会被copied,否则将X将会被覆写。
       n_jobs：默认值为1，指定训练时用的CPU资源。
       '''
       #LinearRegression(fit_intercept=True,normalize=False,copy_X=True,n_jobs=1)
       model=LinearRegression()
       #用训练数据拟合模型
       model.fit(X_train,y_train)
       #输出模型的相关因子和截距
       print(model.coef_,model.intercept_)
       print(model.get_params())#查看模型的参数
       y_hat=model.predict(X_test)
       mse=np.average((y_hat-y_test)**2)
       rmse=np.sqrt(mse)
       print('mse:',mse)
       print('rmse',rmse)
       print("Train R2=:",model.score(X_train,y_train))#返回R^2
       print("Test R2=:",model.score(X_test,y_test))
   
       plt.figure(facecolor='w')
       t=np.arange(len(X_test))
       plt.plot(t,y_test,'r-',label='真实数据',lw=2)
       plt.plot(t,y_hat,'g-',lw=2,label='预测数据')
       plt.legend(loc='upper right')
       plt.title('线性回归预测销量',fontsize=18)
       plt.grid(b=True,ls=':')
       plt.show()

Lasso/Ridge代码

import pandas as pd
from sklearn.linear_model import LinearRegression,Lasso,Ridge
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
from sklearn.model_selection import train_test_split,GridSearchCV
if __name__=='__main__':
    data=pd.read_csv('./Data/Advertising.csv')
    X=data[['TV','Radio','Newspaper']]
    y=data['Sales']
    #将数据集和训练集分开
    X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.3,random_state=1)
    model=Lasso()
    # model=Ridge()
    # model=LinearRegression()
    #超参数的候选值，-3~2选10个值
    alpha=np.logspace(-3,2,10)
    #不用科学计数法
    np.set_printoptions(suppress=True)
    #网格搜索法，为了选取超参数而进行交叉验证
    lasso_model=GridSearchCV(model,param_grid={'alpha':alpha},cv=10)
    #拟合训练数据
    lasso_model.fit(X_train,y_train)
    #打印出模型的所有参数
    print("参数:",lasso_model.get_params())
    #打印出模型的超参数
    print('超参数',lasso_model.best_params_)

    order=y_test.argsort(axis=0)
    y_test=y_test.values[order]
    x_test=X_test.values[order,:]
    y_hat=lasso_model.predict(x_test)
    print(lasso_model.score(x_test,y_test))#返回的是R^2
    mse=np.average((y_hat-y_test)**2)
    rmse=mse*1/2
    print(mse,rmse)

    t=np.arange(len(X_test))
    # 用来显示正常中文标签
    mpl.rcParams['font.sans-serif'] = ['simHei']
    # 用来显示正常负号
    mpl.rcParams['axes.unicode_minus'] = False
    plt.figure(facecolor='w')
    plt.plot(t, y_test, 'r-', label='真实数据', lw=2)
    plt.plot(t, y_hat, 'g-', lw=2, label='预测数据')
    plt.legend(loc='upper right')
    plt.title('线性回归预测销量', fontsize=18)
    plt.grid(b=True, ls=':')
    plt.show()

补充

1. 矩阵对向量求偏导：
   A为$m×n$的矩阵,$\vec x$为$n×1$的列向量，则$A\vec x$为m×1的列向量，记做$\vec y$=$A\vec x$
   思考：$\frac{\partial \vec y}{\partial \vec x}=?$

   $$A\vec x=\begin{bmatrix} a_{11}&a_{12}&\cdots &a_{1n}\\ a_{21}&a_{22}&\cdots &a_{2n}\\ \vdots &\vdots &\ddots &\vdots\\ a_{m1}&a_{m2}&\cdots &a_{mn}\\ \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_m\\ \end{bmatrix} =\begin{bmatrix} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n\\ a_{21}x_2+a_{22}x_2+\cdots+a_{2n}x_n\\ \vdots\\ a_{m1}x_1+a_{m2}x_2+\cdots +a_{mn}x_n\\ \end{bmatrix}$$

   首先考虑矩阵$A$先对$x_1$求偏导得到第一行，以此类推一直得到得到第$m$行,因而，

   $$\frac{\partial \vec y}{\partial \vec x}= \begin{bmatrix} a_{11}&a_{21}&\cdots&a_{n1}\\ a_{12}&a_{22}&\cdots&a_{n2}\\ \vdots &\vdots &\ddots &\vdots\\ a_{1m}&a_{2m}&\cdots&a_{nm}\\ \end{bmatrix}=A^T$$

   结论推广：
           $\frac{\partial A\vec x}{\partial x}=A^T$  $\frac{\partial A\vec x}{\partial x^T}=A$  $\frac{\partial {\vec x}^T A}{\partial x}=A$

2. 标量对向量的导数
   $A$为$n×n$的对称阵($A^T=A$),$\vec x$为$n×1$的列向量,记做$y$=$A\vec x$,则

   $$\frac{\partial y}{\partial \vec x}=\frac{\vec x^T·A·\vec x}{\partial \vec x}=(A^T+A)·\vec x=2A\vec x$$

   推导：

   $$\vec x^TA\vec x=\begin{pmatrix}x_1&x_2&\cdots&x_n\end{pmatrix}\begin{pmatrix}\sum_{j=1}^n a_{1j}x_j&\sum_{j=1}^n a_{2j}x_j&\cdots&\sum_{j=1}^n a_{nj}x_j\end{pmatrix}^T$$ $$=\sum_{i=1}^n((\sum_{j=1}^n a_{ij}x_j )x_i)=\sum_{i=1}^n \sum_{j=1}^n a_{ij}x_ix_j$$

   则得到如下结论，其中一项可以看做是对第i行的求偏导，一项可以看做是对i列的求偏导：

   $$\frac{(\vec x^TA\vec x)}{x_i}=\sum_{j=1}^na_{ij}x_j+\sum_{j=1}^na_{ji}x_j=\sum_{j=1}^n(a_{ij}+a_{ji})x_j$$

3. 标量对方阵求导数：
   $A$为$n×n$的矩阵,$|A|$为$A$的行列式,计算$\frac{\partial |A|}{A}$.
   计算：根据行列式的计算公式，对于$\forall 1\le i\le n$,$|A|=\sum_{j=0}^na_{ij}(-1)^{i+j}M_{ij}$,其中，$M_{ij}$是$a_{ij}$的代数余子式。
   则有：

   $$\frac{\partial |A|}{\partial a_{ij}}=\frac{\partial \sum_{j=0}^na_{ij}(-1)^{i+j}M_{ij} }{\partial a_{ij}}=(-1)^{i+j}M_{ij}=A_{ji}^{\*}$$

   并且$A·A*=|A|·I$,从而:$\frac{\partial |A|}{\partial A}=(A^*)^T=|A|(A^{-1})^T$

   附录

   $x^{(i)}$代表的是样本的特征
   $y^{(i)}$代表的是第$n$个样本的$label$，可离散可连续
   $\theta_i$代表的某个样本中第$i$个数据的参数

   参考资料

   邹博.机器学习
   三种梯度下降方式比较
   L1 L2范数最优化过程-知乎