机器学习（上）-优化基础模型

1. 优化基础模型

回归的基本算法是使用数据集去估计模型的参数，如线性回归模型中的参数w，那么这个数据集我们称为训练数据集，简称训练集。在训练过程中，一般会将损失函数在训练集达到最小值，实际问题往往会让损失函数在训练集最小化为0。

如：在线性回归中，加入非常多的高次项，使模型在训练集的每一个数据点都恰好位于曲线上，那这时候模型在训练集的损失值也就是误差为0。

这个时候也不代表模型可以预测任意情况，因为模型的初衷是为了探索未知/情况复杂的测试数据集（未出现在训练集的未知数据），简称测试集，让其达到最优。（eg:根据股票市场前六个月的数据拟合一个预测模型，目的是预测明天乃至未来的股价变化,不是为了预测以前这六个月越准越好。）

<font color ='red'>**(a) 训练均方误差与测试均方误差：** </font>    

在回归中，最常用的评价指标为均方误差，即：$MSE=\frac{1}{N}\sum _{i=1}^N(y_i-\hat{f}(x_i){)}^2$，其中$\hat{f}(x_i)$是样本$x_i$应用建立的模型$\hat{f}$预测的结果。

均方误差会根据计算数据是用的训练集还是测试集，是训练集上的数据叫训练均方误差，测试集的数据叫测试均方误差。一般模型不会考虑训均方误差，比较关心测试均方误差达到最小。

如何选择一个测试误差最小的模型？

模型建立阶段，不能得到测试数据（比如：我们在模型未上线之前是不能拿到未知且真实的测试数据来验证我们的模型），为了方便，有些会通过训练误差最小化来选择模型，但这个比较片面，且行不通，因为训练误差最小无法让测试误差也最小，有可能会出现测试均方误差很大的情况，这种就是**过拟合**！如图：

(b) 偏差-方差的权衡：

从上图的测试均方误差曲线可以看到：测试均方误差曲线呈现U型曲线，这表明了在测试误差曲线中有两种力量在互相博弈。可以证明：
$$E{\left(y_0-\hat{f}\left(x_0\right)\right)}^2=\mathrm{Var}\left(\hat{f}\left(x_0\right)\right)+{\left[\mathrm{Bias}\left(\hat{f}\left(x_0\right)\right)\right]}^2+\mathrm{Var}(\epsilon )$$
测试均方误差的期望值可以分解为$\hat{f}(x_0)$的方差、$\hat{f}(x_0)$的偏差平方和误差项$ϵ$的方差。

要想模型的测试均方误差达到最小值，也就是同时最小化偏差的平方和方差。但偏差平方和方差本身是非负的，因此测试均方误差的期望不可能会低于误差的方差，因此我们称$\mathrm{Var}(\epsilon )$为建模任务的难度，这个量在我们的任务确定后是无法改变的，也叫做不可约误差。

tips:

模型的方差：用不同的数据集去估计$f$时，估计函数的改变量。

eg：想要建立一个线性回归模型，可以通过输入中国人身高去预测我们的体重。实际肯定无法把全中国13亿人做一次人口普查，拿到他们的身高体重去建立模型。所以，我只能是13亿中抽1000(N)个样本进行建模，然后对这个抽样的过程重复100遍，就会得到100个1000人的样本集。接着使用线性回归模型估计参数就能得到100个线性回归模型。由于样本抽取具有随机性，我们得到的100个模型不可能参数完全一样，那么这100个模型之间的差异就叫做方差。

显然，想得到一个稳定的模型，也就是在不同的样本集估计的模型都不会相差太大，即要求f的方差越小越好。一般来说，模型的复杂度越高，f的方差就会越大。 如加入二次项的模型的方差比线性回归模型的方差要大。

模型的偏差：选择一个简单的模型去估计真实函数所带入的误差。

eg：假如真实的数据X与Y的关系是二次关系，但是我们选择了线性模型进行建模，那由于模型的复杂度引起的这种误差我们称为偏差，它的构成时复杂的。偏差度量了学习算法的期望预测与真实结果的偏离程度，即刻画了学习算法本身的拟合能力。

偏差度量的是单个模型的学习能力，而方差度量的是同一个模型在不同数据集上的稳定性。

“偏差-方差分解”说明：泛化性能是由学习算法的能力、数据的充分性以及学习任务本身的难度所共同决定的。给定学习任务，为了取得好的泛化性能，则需使偏差较小，即能够充分拟合数据，并且使方差较小，即使得数据扰动产生的影响小。

一般而言，增加模型的复杂度，会增加模型的方差，但是会减少模型的偏差，我们要找到一个方差–偏差的权衡，使得测试均方误差最小。

© 特征提取：
要选择一个测试误差达到最小的模型，要对测试误差进行估计，估计的方式有两种：训练误差修正与交叉验证。

训练误差修正：

模型越复杂，训练误差越小，测试误差先减后增。因此，先构造一个特征较多的模型使其过拟合，此时训练误差很小而测试误差很大，那这时我们加入关于特征个数的惩罚。当我们的训练误差随着特征个数的增加而减少时，惩罚项因为特征数量的增加而增大，抑制了训练误差随着特征个数的增加而无休止地减小。具体的数学量如下：
$C_p=\frac{1}{N}(RSS+2d{\hat{\sigma }}^2)$，其中d为模型特征个数，$RSS=\sum _{i=1}^N(y_i-\hat{f}(x_i){)}^2$，${\hat{\sigma }}^2$为模型预测误差的方差的估计值，即残差的方差。
AIC赤池信息量准则：$AIC=\frac{1}{d{\hat{\sigma }}^2}(RSS+2d{\hat{\sigma }}^2)$
BIC贝叶斯信息量准则：$BIC=\frac{1}{n}(RSS+log(n)d{\hat{\sigma }}^2)$

交叉验证：

训练误差修正得到测试误差的估计是间接方法，这种方法的桥梁是训练误差，而交叉验证则是对测试误差的直接估计。
交叉验证比训练误差修正的优势在于：能够给出测试误差的一个直接估计。
在这里只介绍K折交叉验证：我们把训练样本分成K等分，然后用K-1个样本集当做训练集，剩下的一份样本集为验证集去估计由K-1个样本集得到的模型的精度，这个过程重复K次取平均值得到测试误差的一个估计$C{V}_{(K)}=\frac{1}{K}\sum _{i=1}^{K}MS{E}_i$。5折交叉验证如下图：（蓝色的是训练集，黄色的是验证集）

在测试误差能够被合理的估计出来以后，我们做特征选择的目标就是：从p个特征中选择m个特征，使得对应的模型的测试误差的估计最小。对应的方法有：

最优子集选择:

(i) 记不含任何特征的模型为$M_0$，计算这个$M_0$的测试误差。
(ii) 在$M_0$基础上增加一个变量，计算p个模型的RSS，选择RSS最小的模型记作$M_1$，并计算该模型$M_1$的测试误差。
(iii) 再增加变量，计算p-1个模型的RSS，并选择RSS最小的模型记作$M_2$，并计算该模型$M_2$的测试误差。
(iv) 重复以上过程知道拟合的模型有p个特征为止，并选择p+1个模型 { M 0 , M 1 , . . . , M p } \{M_0,M_1,...,M_p \} {M0​,M1​,...,Mp​}中测试误差最小的模型作为最优模型。

向前逐步选择:

最优子集选择虽然在原理上很直观，但是随着数据特征维度p的增加，子集的数量为$2^p$，计算效率非常低下且需要的计算内存也很高，在大数据的背景下显然不适用。因此，我们需要把最优子集选择的运算效率提高，因此向前逐步选择算法的过程如下：
(i) 记不含任何特征的模型为$M_0$，计算这个$M_0$的测试误差。
(ii) 在$M_0$基础上增加一个变量，计算p个模型的RSS，选择RSS最小的模型记作$M_1$，并计算该模型$M_1$的测试误差。
(iii) 在最小的RSS模型下继续增加一个变量，选择RSS最小的模型记作$M_2$，并计算该模型$M_2$的测试误差。
(iv) 以此类推，重复以上过程知道拟合的模型有p个特征为止，并选择p+1个模型 { M 0 , M 1 , . . . , M p } \{M_0,M_1,...,M_p \} {M0​,M1​,...,Mp​}中测试误差最小的模型作为最优模型。

(d) 压缩估计(正则化)：
刚刚是直接对特征自身进行选择，还可以对回归的系数进行约束或者加罚，对p个特征的模型进行拟合，显著降低模型方差，这样也会提高模型的拟合效果。具体来说，就是将回归系数往零的方向压缩，即压缩估计。

**岭回归(L2正则化的例子): **

在线性回归中，我们的损失函数为$J(w)=\sum _{i=1}^N(y_i-w_0-\sum _{j=1}^p{w}_j{x}_{ij}{)}^2$，我们在线性回归的损失函数的基础上添加对系数的约束或者惩罚，即：
$$\begin{gathered} J(w)=\sum _{i=1}^N(y_i-w_0-\sum _{j=1}^p{w}_j{x}_{ij}{)}^2+\lambda \sum _{j=1}^p{w}_j^2,其中，\lambda \ge 0 \\ \hat{w}=(X^{T}X+\lambda I{)}^{-1}{X}^{T}Y \end{gathered}$$
调节参数$\lambda$的大小是影响压缩估计的关键，$\lambda$越大，惩罚的力度越大，系数则越趋近于0，但无法等于0。反之，选择合适的$\lambda$对模型精度来说十分重要。岭回归通过牺牲线性回归的无偏性降低方差，有可能使得模型整体的测试误差较小，提高模型的泛化能力。

Lasso回归(L1正则化的例子)：

岭回归的系数只能趋于0但无法等于0，即无法做特征选择。能否使用压缩估计的思想做到像特征最优子集选择那样提取出重要的特征呢？

只需要对岭回归的优化函数做小小的调整，使用系数向量的L1范数替换岭回归中的L2范数：
$$J(w)=\sum _{i=1}^N(y_i-w_0-\sum _{j=1}^p{w}_j{x}_{ij}{)}^2+\lambda \sum _{j=1}^p|w_j|,其中，\lambda \ge 0$$
为什么Losso能做到特征选择而岭回归却不能做到呢？(如图：左边为lasso，右边为岭回归)

椭圆形曲线为RSS等高线，菱形和圆形区域分别代表了L1和L2约束，Lsaao回归和岭回归都是在约束下的回归，因此最优的参数为椭圆形曲线与菱形和圆形区域相切的点。
但是Lasso回归的约束在每个坐标轴上都有拐角，因此当RSS曲线与坐标轴相交时恰好回归系数中的某一个为0，这样就实现了特征提取。
反观岭回归的约束是一个圆域，没有尖点，因此与RSS曲线相交的地方一般不会出现在坐标轴上，因此无法让某个特征的系数为0，因此无法做到特征提取。

(e) 降维：

目前基于原始特征$x_1,...,x_p$得到的2种方法可以对方差的控制：

①使用原始变量的子集
②将变量系数压缩至零

第③种方法：改变原始特征

将原始的特征空间投影到一个低维的空间实现变量的数量变少。（如：将二维的平面投影至一维空间。）

机器学习领域中的降维：

采用某种映射方法，将原高维空间中的数据点映射到低维度的空间中，将冗余信息以及噪音信息减少，减小误差，提高识别精度，或者寻找数据内部的本质结构特征。

降维的本质是学习一个映射函数 f : x->y

x：原始数据点的表达，目前最多使用向量表达形式。

y：数据点映射后的低维向量表达，通常y的维度小于x的维度（当然提高维度也是可以的）。

f：显式的或隐式的、线性的或非线性的。

降维算法成为了数据预处理的一部分，如PCA。

主成分分析(PCA)：

通过最大投影方差 将原始空间进行重构，即由特征相关重构为无关，即落在某个方向上的点(投影)的方差最大。在进行下一步推导之前，我们先把样本均值和样本协方差矩阵推广至矩阵形式：
样本均值Mean:$\bar{x}=\frac{1}{N}\sum _{i=1}^N{x}_i=\frac{1}{N}{X}^T{1}_N,其中1_N=(1,1,...,1{)}_N^T$
样本协方差矩阵$S^2=\frac{1}{N}\sum _{i=1}^N(x_i-\bar{x})(x_i-\bar{x}{)}^T=\frac{1}{N}{X}^{T}HX,其中，H=I_N-\frac{1}{N}{1}_N{1}_N^T$

最大投影方差的步骤：

(i) 中心化：$x_i-\bar{x}$
(ii) 计算每个点$x_1,...,x_N$至${\vec{u}}_1$方向上的投影：$(x_i-\bar{x}){\vec{u}}_1,||{\vec{u}}_1||=1$
(iii) 计算投影方差：$J=\frac{1}{N}\sum _{i=1}^N[(x_i-\bar{x}{)}^T{\vec{u}}_1{]}^2,||{\vec{u}}_1||=1$
(iv) 最大化投影方差求${\vec{u}}_1$：
$$\begin{gathered} {\bar{u}}_1=argma{x}_{u_1}\frac{1}{N}\sum _{i=1}^N[(x_i-\bar{x}{)}^T{\vec{u}}_1{]}^2 \\ s.t.{\vec{u}}_1^T{\vec{u}}_1=1({\vec{u}}_1往后不带向量符号) \end{gathered}$$
得到：
$$\begin{gathered} J=\frac{1}{N}\sum _{i=1}^N[(x_i-\bar{x}{)}^T{\vec{u}}_1{]}^2=\frac{1}{N}\sum _{i=1}^N[u_1^T(x_i-\bar{x})(x_i-\bar{x}{)}^T{u}_1] \\ =u_1^T[\frac{1}{N}\sum _{i=1}^N(x_i-\bar{x})(x_i-\bar{x}{)}^T]u_1=u_1^T{S}^2{u}_1 \end{gathered}$$
即：
$$\begin{gathered} {\hat{u}}_1=argma{x}_{u_1}{u}_1^T{S}^2{u}_1,s.t.u_1^T{u}_1=1 \\ L(u_1,\lambda )=u_1^T{S}^2{u}_1+\lambda (1-u_1^T{u}_1) \\ \frac{\partial L}{\partial u_1}=2S^2{u}_1-2\lambda u_1=0 \\ 即：S^2{u}_1=\lambda u_1 \end{gathered}$$
可以看到：$\lambda$为$S^2$的特征值，$u_1$为$S^2$的特征向量。因此我们只需要对中心化后的协方差矩阵进行特征值分解，得到的特征向量即为投影方向。如果需要进行降维，那么只需要取p的前M个特征向量即可。

# 引入相关科学计算包
import numpy as np
import pandas as pd
# import matplotlib.pyplot as plt 
# %matplotlib inline 
# plt.style.use("ggplot")  
# import seaborn as sns

from sklearn import datasets
boston = datasets.load_boston()     # 返回一个类似于字典的类
X = boston.data
y = boston.target
features = boston.feature_names
boston_data = pd.DataFrame(X,columns=features)
boston_data["Price"] = y
boston_data.head()

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT	Price
0	0.00632	18.0	2.31	0.0	0.538	6.575	65.2	4.0900	1.0	296.0	15.3	396.90	4.98	24.0
1	0.02731	0.0	7.07	0.0	0.469	6.421	78.9	4.9671	2.0	242.0	17.8	396.90	9.14	21.6
2	0.02729	0.0	7.07	0.0	0.469	7.185	61.1	4.9671	2.0	242.0	17.8	392.83	4.03	34.7
3	0.03237	0.0	2.18	0.0	0.458	6.998	45.8	6.0622	3.0	222.0	18.7	394.63	2.94	33.4
4	0.06905	0.0	2.18	0.0	0.458	7.147	54.2	6.0622	3.0	222.0	18.7	396.90	5.33	36.2

特征提取的实例：向前逐步回归
案例来源：https://blog.csdn.net/weixin_44835596/article/details/89763300
根据AIC准则定义向前逐步回归进行变量筛选

#定义向前逐步回归函数
def forward_select(data,target):
    variate=set(data.columns)  #将字段名转换成字典类型
    variate.remove(target)  #去掉因变量的字段名
    selected=[]
    current_score,best_new_score=float('inf'),float('inf')  #目前的分数和最好分数初始值都为无穷大（因为AIC越小越好）
    #循环筛选变量
    while variate:
        aic_with_variate=[]
        for candidate in variate:  #逐个遍历自变量
            formula="{}~{}".format(target,"+".join(selected+[candidate]))  #将自变量名连接起来
            aic=ols(formula=formula,data=data).fit().aic  #利用ols训练模型得出aic值
            aic_with_variate.append((aic,candidate))  #将第每一次的aic值放进空列表
        aic_with_variate.sort(reverse=True)  #降序排序aic值
        best_new_score,best_candidate=aic_with_variate.pop()  #最好的aic值等于删除列表的最后一个值，以及最好的自变量等于列表最后一个自变量
        if current_score>best_new_score:  #如果目前的aic值大于最好的aic值
            variate.remove(best_candidate)  #移除加进来的变量名，即第二次循环时，不考虑此自变量了
            selected.append(best_candidate)  #将此自变量作为加进模型中的自变量
            current_score=best_new_score  #最新的分数等于最好的分数
            print("aic is {},continuing!".format(current_score))  #输出最小的aic值
        else:
            print("for selection over!")
            break
    formula="{}~{}".format(target,"+".join(selected))  #最终的模型式子
    print("final formula is {}".format(formula))
    model=ols(formula=formula,data=data).fit()
    return(model)

import statsmodels as sm #最小二乘法
from statsmodels.formula.api import ols #加载ols模型
forward_select(data=boston_data,target="Price")

aic is 3286.974956900157,continuing!
aic is 3171.5423142992013,continuing!
aic is 3114.0972674193326,continuing!
aic is 3097.359044862759,continuing!
aic is 3069.438633167217,continuing!
aic is 3057.9390497191152,continuing!
aic is 3048.438382711162,continuing!
aic is 3042.274993098419,continuing!
aic is 3040.166722868547,continuing!
aic is 3032.111008343286,continuing!
aic is 3021.7984119907965,continuing!
for selection over!
final formula is Price~LSTAT+RM+PTRATIO+DIS+NOX+CHAS+B+ZN+CRIM+RAD+TA
<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1b233f2fb00>

此处调用import statsmodels.api as sm报错：ImportError: cannot import name ‘datetools’

解决方法：更新statsmodels库，一定要关掉程序，然后重新启动，再运行

pip install -U --pre statsmodels

如果还报错

解决方法：高版本的在import的时候不需要加.api

lm=ols("Price~LSTAT+RM+PTRATIO+DIS+NOX+CHAS+B+ZN+CRIM+RAD+TAX",data=boston_data).fit()
lm.summary()

OLS Regression Results

Dep. Variable:	Price	R-squared:	0.741
Model:	OLS	Adj. R-squared:	0.735
Method:	Least Squares	F-statistic:	128.2
Date:	Thu, 18 Mar 2021	Prob (F-statistic):	5.74e-137
Time:	11:15:40	Log-Likelihood:	-1498.9
No. Observations:	506	AIC:	3022.
Df Residuals:	494	BIC:	3073.
Df Model:	11
Covariance Type:	nonrobust

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	36.3694	5.069	7.176	0.000	26.411	46.328
LSTAT	-0.5232	0.047	-11.037	0.000	-0.616	-0.430
RM	3.7966	0.406	9.343	0.000	2.998	4.595
PTRATIO	-0.9471	0.129	-7.337	0.000	-1.201	-0.693
DIS	-1.4934	0.186	-8.039	0.000	-1.858	-1.128
NOX	-17.3956	3.536	-4.920	0.000	-24.343	-10.448
CHAS	2.7212	0.854	3.185	0.002	1.043	4.400
B	0.0094	0.003	3.508	0.000	0.004	0.015
ZN	0.0458	0.014	3.387	0.001	0.019	0.072
CRIM	-0.1076	0.033	-3.296	0.001	-0.172	-0.043
RAD	0.2991	0.063	4.719	0.000	0.175	0.424
TAX	-0.0118	0.003	-3.488	0.001	-0.018	-0.005

Omnibus:	178.444	Durbin-Watson:	1.078
Prob(Omnibus):	0.000	Jarque-Bera (JB):	786.944
Skew:	1.524	Prob(JB):	1.31e-171
Kurtosis:	8.295	Cond. No.	1.47e+04

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.47e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

岭回归实例分享：
sklearn.linear_model.ridge_regression(X, y, alpha, *, sample_weight=None, solver=‘auto’, max_iter=None, tol=0.001, verbose=0, random_state=None, return_n_iter=False, return_intercept=False, check_input=True)
https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ridge_regression.html?highlight=rid#sklearn.linear_model.ridge_regression

• 参数：
  alpha：较大的值表示更强的正则化。浮点数
  sample_weight：样本权重，默认无。
  solver：求解方法，{‘auto’, ‘svd’, ‘cholesky’, ‘lsqr’, ‘sparse_cg’, ‘sag’, ‘saga’}, 默认=’auto’。“ svd”使用X的奇异值分解来计算Ridge系数。'cholesky’使用标准的scipy.linalg.solve函数通过dot（XT，X）的Cholesky分解获得封闭形式的解。'sparse_cg’使用scipy.sparse.linalg.cg中的共轭梯度求解器。作为一种迭代算法，对于大规模数据（可能设置tol和max_iter），此求解器比“ Cholesky”更合适。 lsqr”使用专用的正则化最小二乘例程scipy.sparse.linalg.lsqr。它是最快的，并且使用迭代过程。“ sag”使用随机平均梯度下降，“ saga”使用其改进的无偏版本SAGA。两种方法都使用迭代过程，并且当n_samples和n_features都很大时，通常比其他求解器更快。请注意，只有在比例大致相同的要素上才能确保“ sag”和“ saga”快速收敛。您可以使用sklearn.preprocessing中的缩放器对数据进行预处理。最后五个求解器均支持密集和稀疏数据。但是，当fit_intercept为True时，仅’sag’和’sparse_cg’支持稀疏输入。

from sklearn import linear_model
reg_rid = linear_model.Ridge(alpha=.5)
reg_rid.fit(X,y)
reg_rid.score(X,y)

0.739919756983817

Lasso实例分享：
class sklearn.linear_model.Lasso(alpha=1.0, *, fit_intercept=True, normalize=False, precompute=False, copy_X=True, max_iter=1000, tol=0.0001, warm_start=False, positive=False, random_state=None, selection=‘cyclic’)
https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html?highlight=lasso#sklearn.linear_model.Lasso

• 参数：
  alpha：正则化强度，1.0代表标准最小二乘。
  fit_intercept：是否计算模型截距。默认true。
  normalize：是否标准化，默认false。
  positive：是否强制系数为正，默认false。

from sklearn import linear_model
reg_lasso = linear_model.Lasso(alpha = 0.5)
reg_lasso.fit(X,y)
reg_lasso.score(X,y)

0.7139506197036717