机器学习基础-sklearn全流程实战

本文转载自DateWhale集成学习，增加了一些个人学习过程的注释，有需要的同学可以参考原文
原文链接
视频地址

1.导论

什么是机器学习？机器学习的一个重要的目标就是利用数学模型来理解数据，发现数据中的规律，用作数据的分析和预测。数据通常由一组向量组成，这组向量中的每个向量都是一个样本，我们用$x_i$来表示一个样本，其中$i=1,2,3,...,N$,共N个样本，每个样本$x_i=(x_{i1},x_{i2},...,x_{ip},y_i)$共p+1个维度，前p个维度的每个维度我们称为一个特征，最后一个维度$y_i$我们称为因变量(响应变量)。特征用来描述影响因变量的因素，如：我们要探寻身高是否会影响体重的关系的时候，身高就是一个特征，体重就是一个因变量。通常在一个数据表dataframe里面，一行表示一个样本$x_i$，一列表示一个特征。
根据数据是否有因变量，机器学习的任务可分为：有监督学习和无监督学习。

• 有监督学习：给定某些特征去估计因变量，即因变量存在的时候，我们称这个机器学习任务为有监督学习。如：我们使用房间面积，房屋所在地区，环境等级等因素去预测某个地区的房价。
• 无监督学习：给定某些特征但不给定因变量，建模的目的是学习数据本身的结构和关系。如：我们给定某电商用户的基本信息和消费记录，通过观察数据中的哪些类型的用户彼此间的行为和属性类似，形成一个客群。注意，我们本身并不知道哪个用户属于哪个客群，即没有给定因变量。
  

根据因变量的是否连续，有监督学习又分为回归和分类：

• 回归：因变量是连续型变量，如：房价，体重等。
• 分类：因变量是离散型变量，如：是否患癌症，西瓜是好瓜还是坏瓜等。

为了更好地叙述后面的内容，我们对数据的形式作出如下约定：
第i个样本：$x_i=(x_{i1},x_{i2},...,x_{ip},y_i{)}^T,i=1,2,...,N$
因变量$y=(y_1,y_2,...,y_N{)}^T$
第k个特征:$x^{(k)}=(x_{1k},x_{2k},...,x_{Nk}{)}^T$
特征矩阵$X=(x_1,x_2,...,x_N{)}^T$

在学习机器学习中，我们经常使用scikit-learn简称sklearn工具库来探索机器学习项目，下面我们开始使用sklearn来演示这几个具体的概念：

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

1.1 回归

首先，我们先来看看有监督学习中回归的例子，我们使用sklearn内置数据集Boston房价数据集。sklearn中所有内置数据集都封装在datasets对象内：
返回的对象有：

• data:特征X的矩阵(ndarray)
• target:因变量的向量(ndarray)
• feature_names:特征名称(ndarray)

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()

sns.scatterplot(boston_data['NOX'],boston_data['Price'],color="r",alpha=0.6)
plt.title("Price~NOX")
plt.show()

可以从数据中看出，因变量Price是一个连续型变量，所以这是一个回归的案例。

1.2 分类

我们来看看大名鼎鼎的iris数据集：

from sklearn import datasets
iris = datasets.load_iris()
X = iris.data
y = iris.target
features = iris.feature_names
iris_data = pd.DataFrame(X,columns=features)
iris_data['target'] = y
iris_data.info()

#可视化特征
marker = ['s','x','o']
for index,c in enumerate(np.unique(y)):  #将一个可遍历的数据对象(如列表、元组或字符串)组合为一个索引序列，同时列出数据和数据下标
    #c=0\1\2，y==c查找与target相同的行，0、1、2各50行
    plt.scatter(x=iris_data.loc[y==c,"sepal length (cm)"],y=iris_data.loc[y==c,"sepal width (cm)"],alpha=0.8,label=c,marker=marker[c])
plt.xlabel("sepal length (cm)")
plt.ylabel("sepal width (cm)")
plt.legend()
plt.show()

我们可以看到：每种不同的颜色和点的样式为一种类型的鸢尾花，数据集有三种不同类型的鸢尾花。因此因变量是一个类别变量，因此通过特征预测鸢尾花类别的问题是一个分类问题。

各个特征的相关解释：

• sepal length (cm)：花萼长度(厘米)
• sepal width (cm)：花萼宽度(厘米)
• petal length (cm)：花瓣长度(厘米)
• petal width (cm)：花瓣宽度(厘米)

1.3 无监督学习

我们可以使用sklearn生成符合自身需求的数据集，下面我们用其中几个函数例子来生成无因变量的数据集： 更多数据集点击这里

# 生成月牙型非凸集
from sklearn import datasets
x, y = datasets.make_moons(n_samples=2000, shuffle=True,
                  noise=0.05, random_state=None)
for index,c in enumerate(np.unique(y)):
    plt.scatter(x[y==c,0],x[y==c,1],s=7)
plt.show()

# 生成符合正态分布的聚类数据
from sklearn import datasets
x, y = datasets.make_blobs(n_samples=5000, n_features=2, centers=3)
for index,c in enumerate(np.unique(y)):
    plt.scatter(x[y==c, 0], x[y==c, 1],s=7)
plt.show()

在本文中，我们要重点介绍回归和分类的集成学习的问题，因此我们在接下来的章节中不再介绍关于无监督学习的具体算法，后面的内容仅仅涉及回归和分类问题。

2. 使用sklearn构建完整的机器学习项目流程

一般来说，一个完整的机器学习项目分为以下步骤：

• 明确项目任务：回归/分类
• 收集数据集并选择合适的特征。
• 选择度量模型性能的指标。
• 选择具体的模型并进行训练以优化模型。
• 评估模型的性能并调参。

2.1 明确项目任务：回归/分类：

本次任务是预测Boston的房价数据集，由于房价是连续型变量，所以此次任务为回归。

2.2 收集数据集并选择合适的特征：

在数据集上我们使用我们比较熟悉的Boston房价数据集，原因是：

• 第一个，我们通过这些简单的数据集快速让我们上手sklearn，以及掌握sklearn的相关操作。
• 第二个，我们用简单的数据集能更加清晰地介绍机器学习的相关模型，避免在处理数据上花费较大的精力。

import pandas as pd
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：规划地段超过25,000平方英尺的住宅用地比例
• INDUS：城镇非零售商业用地比例
• CHAS：是否在查尔斯河边(=1是)
• NOX：一氧化氮浓度(/千万分之一)
• RM：每个住宅的平均房间数
• AGE：1940年以前建造的自住房屋的比例
• DIS：到波士顿五个就业中心的加权距离
• RAD：放射状公路的可达性指数
• TAX：全部价值的房产税率(每1万美元)
• PTRATIO：按城镇分配的学生与教师比例
• B：1000(Bk - 0.63)^2其中Bk是每个城镇的黑人比例
• LSTAT：较低地位人口
• Price：房价

2.3 选择度量模型性能的指标：

注意，分类和回归模型尤其专属的度量指标。sklearn回归度量-英文 sklearn回归度量-中文

• MSE均方误差：$\text{MSE}(y,\hat{y})=\frac{1}{n_{\text{samples}}}{\sum }_{i=0}^{n_{\text{samples}}-1}(y_i-{\hat{y}}_i{)}^2.$
• MAE平均绝对误差:$\text{MAE}(y,\hat{y})=\frac{1}{n_{\text{samples}}}{\sum }_{i=0}^{n_{\text{samples}}-1}\left|y_i-{\hat{y}}_i\right|$
• $R^2$决定系数：$R^2(y,\hat{y})=1-\frac{{\sum }_{i=1}^n(y_i-{\hat{y}}_i{)}^2}{{\sum }_{i=1}^n(y_i-\bar{y}{)}^2}$ 在计量经济学和线性回归中较多使用
• 解释方差得分: e x p l a i n e d _ v a r i a n c e ( y , y ^ ) = 1 − V a r { y − y ^ } V a r { y } explained\_{}variance(y, \hat{y}) = 1 - \frac{Var\{ y - \hat{y}\}}{Var\{y\}} explained_variance(y,y^​)=1−Var{y}Var{y−y^​}​

在这个案例中，我们使用MSE均方误差为模型的性能度量指标。

2.4 选择具体的模型并进行训练

2.4.1 线性回归模型：

        回归这个概念是19世纪80年代由英国统计学家郎西斯.高尔顿在研究父子身高关系提出来的，他发现：在同一族群中，子代的平均身高介于父代的身高以及族群的平均身高之间。具体而言，高个子父亲的儿子的身高有低于其父亲身高的趋势，而矮个子父亲的儿子身高则有高于父亲的身高的趋势。也就是说，子代的身高有向族群平均身高"平均"的趋势，这就是统计学上"回归"的最初含义。

        回归分析是一种预测性的建模技术，它研究的是因变量（目标）和自变量（特征）之间的关系。这种技术通常用于预测分析，时间序列模型以及发现变量之间的因果关系。通常使用曲线/线来拟合数据点，目标是使曲线到数据点的距离差异最小。而线性回归就是回归问题中的一种，线性回归假设目标值与特征之间线性相关，即满足一个多元一次方程。通过构建损失函数，来求解损失函数最小时的参数w ：
假设：数据集$D=\{(x_1,y_1),...,(x_N,y_N)\}$，$x_i\in R^p,y_i\in R,i=1,2,...,N$，$X=(x_1,x_2,...,x_N{)}^T,Y=(y_1,y_2,...,y_N{)}^T$
假设$X$和$Y$之间存在线性关系，模型的具体形式为$\hat{y}=f(w)=w^{T}x$ 这里假设w和x均为列向量，则y为一个数

图中所示为三维空间，即因变量Y收到两个自变量的影响X1、X2的影响，拟合的线变为平面，即要求平面上下的点到平面的均方误差最小。

可以从下面三个视角去解读如何使得均方误差最小。

2.4.1.1 最小二乘估计：

我们需要衡量真实值$y_i$与线性回归模型的预测值$w^T{x}_i$之间的差距，在这里我们和使用二范数的平方和L(w)来描述这种差距（这里的误差不能用绝对值来表示，因为|x|在x=0不可导），即：
$$\begin{gathered} L(w)=\sum _{i=1}^N||w^T{x}_i-y_i|{|}_2^2=\sum _{i=1}^N(w^T{x}_i-y_i{)}^2=(w^T{X}^T-Y^T)(w^T{X}^T-Y^T{)}^T=w^T{X}^{T}Xw-2w^T{X}^{T}Y+Y{Y}^T \\ 因此，我们需要找到使得L(w)最小时对应的参数w，即： \\ \hat{w}=argminL(w) \\ 为了达到求解最小化L(w)问题，我们应用高等数学的知识，使用求导来解决这个问题： \\ \frac{\partial L(w)}{\partial w}=2X^{T}Xw-2X^{T}Y=0,因此： \\ \hat{w}=(X^{T}X{)}^{-1}{X}^{T}Y \end{gathered}$$

2.4.1.2 几何解释：

在线性代数中，我们知道两个向量a和b相互垂直可以得出：$<a,b>=a.b=a^{T}b=0$,而平面X的法向量为$Y-Xw$，与平面X互相垂直，因此根据两向量a和b相互垂直：$X^T(Y-Xw)=0$，即：$w=(X^{T}X{)}^{-1}{X}^{T}Y$

X1、X2为两个特征，XW为预测值，也肯定位于X平面上，使得真实值Y与预测值$XW$距离最小，即平面X的法向量为Y-Xw

2.4.1.3 概率视角：

假设噪声$ϵ\backsim N(0,{\sigma }^2),y=f(w)+ϵ=w^{T}x+ϵ$，因此：$y|x_i,wN(w^{T}x,{\sigma }^2)$
我们使用极大似然估计MLE对参数w进行估计：
$$\begin{gathered} L(w)=logP(Y|X;w)=log\prod _{i=1}^{N}P(y_i|x_i;w)=\sum _{i=1}^{N}logP(y_i|x_i;w) \\ =\sum _{i=1}^{N}log(\frac{1}{\sqrt{2\pi \sigma }}exp(-\frac{(y_i-w^T{x}_i{)}^2}{2{\sigma }^2}))=\sum _{i=1}^N[log(\frac{1}{\sqrt{2\pi }\sigma })-\frac{1}{2{\sigma }^2}(y_i-w^T{x}_i{)}^2]，因为log(\frac{1}{\sqrt{2\pi }\sigma })为常数 \\ argma{x}_{w}L(w)=argmi{n}_w[l(w)=\sum _{i=1}^N(y_i-w^T{x}_i{)}^2]这个式子正好等于最小二乘法的结果 \\ 因此：线性回归的最小二乘估计<==>噪声ϵ\backsim N(0,{\sigma }^2)的极大似然估计 \end{gathered}$$

下面，我们使用sklearn的线性回归实例来演示：

from sklearn import linear_model      # 引入线性回归方法
lin_reg = linear_model.LinearRegression()       # 创建线性回归的类
lin_reg.fit(X,y)        # 输入特征X和因变量y进行训练
print("模型系数：",lin_reg.coef_)             # 输出模型的系数
print("模型得分：",lin_reg.score(X,y))    # 输出模型的决定系数R^2

模型系数： [-1.08011358e-01  4.64204584e-02  2.05586264e-02  2.68673382e+00
 -1.77666112e+01  3.80986521e+00  6.92224640e-04 -1.47556685e+00
  3.06049479e-01 -1.23345939e-02 -9.52747232e-01  9.31168327e-03
 -5.24758378e-01]
模型得分： 0.7406426641094095

2.4.2 线性回归的推广

在线性回归中，我们假设因变量与特征之间的关系是线性关系，这样的假设使得模型很简单，但是缺点也是显然的，那就是当数据存在非线性关系时，我们使用线性回归模型进行预测会导致预测性能极其低下，因为模型的形式本身是线性的，无法表达数据中的非线性关系。我们一个很自然的想法就是去推广线性回归模型，使得推广后的模型更能表达非线性的关系。

2.4.2.1 多项式回归：

为了体现因变量和特征的非线性关系，一个很自然而然的想法就是将标准的线性回归模型：
$$y_i=w_0+w_1{x}_i+{ϵ}_i$$

换成一个多项式函数：
$$y_i=w_0+w_1{x}_i+w_2{x}_i^2+...+w_d{x}_i^d+ϵ$$
对于多项式的阶数d不能取过大，一般不大于3或者4，因为d越大，多项式曲线就会越光滑，在X的边界处有异常的波动。（图中的边界处的4阶多项式拟合曲线的置信区间(虚线表示置信区间)明显增大，预测效果的稳定性下降。）

左侧的图是阶数比较低的情况，右侧的图是阶数较高的情况。蓝色虚线是置信区间，越大表示越不稳定 。多项式回归实例介绍：
sklearn.preprocessing.PolynomialFeatures(degree=2, *, interaction_only=False, include_bias=True, order=‘C’):

• 参数：
  degree：特征转换的阶数。
  interaction_onlyboolean：是否只包含交互项，默认False 。
  include_bias：是否包含截距项，默认True。
  order：str in {‘C’, ‘F’}, default ‘C’，输出数组的顺序。

from sklearn.preprocessing import PolynomialFeatures
X_arr = np.arange(6).reshape(3, 2)
print("原始X为：\n",X_arr)

poly = PolynomialFeatures(2)
print("2次转化X：\n",poly.fit_transform(X_arr))

poly = PolynomialFeatures(interaction_only=True)
print("2次转化X：\n",poly.fit_transform(X_arr))

2.4.2.2 广义可加模型(GAM)：

广义可加模型GAM实际上是线性模型推广至非线性模型的一个框架，在这个框架中，每一个变量都用一个非线性函数来代替，但是模型本身保持整体可加性。GAM模型不仅仅可以用在线性回归的推广，还可以将线性分类模型进行推广。具体的推广形式是：
标准的线性回归模型：
$$y_i=w_0+w_1{x}_{i1}+...+w_p{x}_{ip}+{ϵ}_i$$
GAM模型框架：
$$y_i=w_0+\sum _{j=1}^p{f}_j(x_{ij})+{ϵ}_i$$
GAM模型的优点与不足：
- 优点：简单容易操作，能够很自然地推广线性回归模型至非线性模型，使得模型的预测精度有所上升；由于模型本身是可加的，因此GAM还是能像线性回归模型一样把其他因素控制不变的情况下单独对某个变量进行推断，极大地保留了线性回归的易于推断的性质。
-

• 推断是指，在其他变量控制不变的情况下分析某个变量对因变量的影响，例如
  $$y_i=w_0+w_1{x}_{i1}$$
  当xi变化1且其他变量不变时，y变化w1，因此就推断出w1的含义。模型越复杂，越难推断。参考《回归分析 谢宇》

• 缺点：GAM模型会经常忽略一些有意义的交互作用，比如某两个特征共同影响因变量，不过GAM还是能像线性回归一样加入交互项$x^{(i)}\times x^{(j)}$的形式进行建模；但是GAM模型本质上还是一个可加模型，如果我们能摆脱可加性模型形式，可能还会提升模型预测精度，详情请看后面的算法。

GAM模型实例介绍：
安装pygam：pip install pygam （sklearn没有集成GAM）

from pygam import LinearGAM
gam = LinearGAM().fit(boston_data[boston.feature_names], y)
gam.summary()

看懂下图需要一定的统计学功底，主要看右侧红色矩阵里的，***越多，表示该变量对因变量的影响越显著。统计学入门《女士品茶》

2.4.3 回归树：

基于树的回归方法主要是依据分层和分割的方式将特征空间划分为一系列简单的区域。对某个给定的待预测的自变量，用他所属区域中训练集的平均数或者众数对其进行预测。由于划分特征空间的分裂规则可以用树的形式进行概括，因此这类方法称为决策树方法。决策树由结点(node)和有向边(diredcted edge)组成。结点有两种类型：内部结点(internal node)和叶结点(leaf node)。内部结点表示一个特征或属性，叶结点表示一个类别或者某个值。区域$R_1,R_2$等称为叶节点，将特征空间分开的点为内部节点。

下图右侧为回归树的模型（例如，对X1的判断，必须贯穿X2，中间不能中断）

建立回归树的过程大致可以分为以下两步：

• 将自变量的特征空间(即$x^{(1)},x^{(2)},x^{(3)},...,x^{(p)}$)的可能取值构成的集合分割成J个互不重叠的区域$R_1,R_2,...,R_j$。

  • 对落入区域$R_j$的每个观测值作相同的预测，预测值等于$R_j$上训练集的因变量的简单算术平均。
    具体来说，就是：
    a. 选择最优切分特征j以及该特征上的最优点s：
    遍历特征j以及固定j后遍历切分点s，选择使得下式最小的(j,s) $mi{n}_{j,s}[mi{n}_{c_1}\sum _{x_i\in R_1(j,s)}(y_i-c_1{)}^2+mi{n}_{c_2}\sum _{x_i\in R_2(j,s)}(y_i-c_2{)}^2]$
    b. 按照(j,s)分裂特征空间： R 1 ( j , s ) = { x ∣ x j ≤ s } 和 R 2 ( j , s ) = { x ∣ x j > s } , c ^ m = 1 N m ∑ x ∈ R m ( j , s ) y i ,    m = 1 , 2 R_1(j,s) = \{x|x^{j} \le s \}和R_2(j,s) = \{x|x^{j} > s \},\hat{c}_m = \frac{1}{N_m}\sum\limits_{x \in R_m(j,s)}y_i,\;m=1,2 R1​(j,s)={x∣xj≤s}和R2​(j,s)={x∣xj>s},c^m​=Nm​1​x∈Rm​(j,s)∑​yi​,m=1,2
    c. 继续调用步骤1，2直到满足停止条件，就是每个区域的样本数小于等于5。
    d. 将特征空间划分为J个不同的区域，生成回归树：$f(x)=\sum _{m=1}^J{\hat{c}}_{m}I(x\in R_m)$
    如以下生成的关于运动员在棒球大联盟数据的回归树：
    

    回归树与线性模型的比较：
    线性模型的模型形式与树模型的模型形式有着本质的区别，具体而言，线性回归对模型形式做了如下假定：$f(x)=w_0+\sum _{j=1}^p{w}_j{x}^{(j)}$，而回归树则是$f(x)=\sum _{m=1}^J{\hat{c}}_{m}I(x\in R_m)$。那问题来了，哪种模型更优呢？这个要视具体情况而言，如果特征变量与因变量的关系能很好的用线性关系来表达，那么线性回归通常有着不错的预测效果，拟合效果则优于不能揭示线性结构的回归树。反之，如果特征变量与因变量的关系呈现高度复杂的非线性，那么树方法比传统方法更优。
    

    树模型的优缺点：

    • 树模型的解释性强，在解释性方面可能比线性回归还要方便。
    • 树模型更接近人的决策方式。
    • 树模型可以用图来表示，非专业人士也可以轻松解读。
    • 树模型可以直接做定性的特征而不需要像线性回归一样哑元化。（参照视频P8）
    • 树模型能很好处理缺失值和异常值，对异常值不敏感，但是这个对线性模型来说却是致命的。
    • 树模型的预测准确性一般无法达到其他回归模型的水平，但是改进的方法很多。（很容易过拟合）

sklearn使用回归树的实例：官网
sklearn.tree.DecisionTreeRegressor(*, criterion=‘mse’, splitter=‘best’, max_depth=None, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features=None, random_state=None, max_leaf_nodes=None, min_impurity_decrease=0.0, min_impurity_split=None, presort=‘deprecated’, ccp_alpha=0.0）

• 参数：(列举几个重要的，常用的，详情请看上面的官网)
  criterion：{“ mse”，“ friedman_mse”，“ mae”}，默认=“ mse”。衡量分割标准的函数 。
  splitter：{“best”, “random”}, default=”best”。分割方式。
  max_depth：树的最大深度。
  min_samples_split：拆分内部节点所需的最少样本数，默认是2。
  min_samples_leaf：在叶节点处需要的最小样本数。默认是1。
  min_weight_fraction_leaf：在所有叶节点处（所有输入样本）的权重总和中的最小加权分数。如果未提供sample_weight，则样本的权重相等。默认是0。

from sklearn.tree import DecisionTreeRegressor    
reg_tree = DecisionTreeRegressor(criterion = "mse",min_samples_leaf = 5)
reg_tree.fit(X,y)
reg_tree.score(X,y)
#0.9376307599929274

2.4.4 支持向量机回归(SVR)

在介绍支持向量回归SVR之前，我们先来了解下约束优化的相关知识：

• 约束优化问题§：

$$\begin{gathered} minf(x) \\ s.t.g_i(x)\le 0,i=1,2,...,m \\ {h}_j(x)=0,j=1,2,...,l \end{gathered}$$

我们假设$x^{\ast }$为满足以上条件的局部最优解，$p^{\ast }=f(x^{\ast })$，我们的目的就是要找到$x^{\ast }$与$p^{\ast }$，满足不等式和等式约束的x集合成为可行域，记作S。

• KKT条件(最优解的一阶必要条件)

因为KKT条件是最优化的相关内容，在本次开源学习中并不是重点，因此在这里我用一个更加简单的例子说明KKT条件，严格的证明请参见凸优化相关书籍。

在这个例子中，我们考虑：($x^{\ast }$为我们的最优解)
$$\begin{gathered} minf(x) \\ s.t.g_1(x)\le 0,x\in R^n \\ {g}_2(x)\le 0 \\ {g}_3(x)\le 0 \end{gathered}$$

我们可以看到：$-\nabla f(x^{\ast })$可以由$\nabla g_1(x^{\ast })$与$\nabla g_2(x^{\ast })$线性表出，因此有：$-\nabla f(x^{\ast })={\lambda }_1\nabla g_1(x^{\ast })+{\lambda }_2\nabla g_2(x^{\ast })$，其中${\lambda }_1,{\lambda }_2\ge 0$，即：
$$\nabla f(x^{\ast })+{\lambda }_1\nabla g_1(x^{\ast })+{\lambda }_2\nabla g_2(x^{\ast })=0,其中{\lambda }_1,{\lambda }_2\ge 0$$
我们把没有起作用的约束$g_3(x)$也放到式子里面去，目的也就是为了书写方便，即要求：
$$\nabla f(x^{\ast })+{\lambda }_1\nabla g_1(x^{\ast })+{\lambda }_2\nabla g_2(x^{\ast })+{\lambda }_3\nabla g_3(x^{\ast })=0,其中{\lambda }_1,{\lambda }_2\ge 0,{\lambda }_3=0$$
由于点$x^{\ast }$位于方程$g_1(x)=0$与$g_2(x)=0$上，因此：${\lambda }_1{g}_1(x^{\ast })=0,{\lambda }_2{g}_2(x^{\ast })=0,{\lambda }_3{g}_3(x^{\ast })=0$

因此，KKT条件就是：假设$x^{\ast }$为最优化问题§的局部最优解，且$x^{\ast }$ 在某个适当的条件下 ,有：
$$\begin{gathered} \nabla f(x^{\ast })+\sum _{i=1}^m{\lambda }_i\nabla g(x^{\ast })+\sum _{j=1}^l{\mu }_j\nabla h_j(x^{\ast })=0(对偶条件) \\ {\lambda }_i\ge 0,i=1,2,...,m(对偶条件) \\ {g}_i(x^{\ast })\le 0(原问题条件) \\ {h}_j(x^{\ast })=0(原问题条件) \\ {\lambda }_{i}g(x^{\ast })=0(互补松弛定理) \end{gathered}$$

• 对偶理论：
  为什么要引入对偶问题呢？是因为原问题与对偶问题就像是一个问题两个角度去看，如利润最大与成本最低等。有时侯原问题上难以解决，但是在对偶问题上就会变得很简单。再者，任何一个原问题在变成对偶问题后都会变成一个凸优化的问题，这点我们后面会有介绍。下面我们来引入对偶问题：
  首先，我们的原问题§是：
  $$\begin{gathered} minf(x) \\ s.t.g_i(x)\le 0,i=1,2,...,m \\ {h}_j(x)=0,j=1,2,...,l \end{gathered}$$
  引入拉格朗日函数：$L(x,\lambda ,\mu )=f(x)+\sum _{i=1}^m{\lambda }_i{g}_i(x)+\sum _{j=1}^l{\mu }_j{h}_j(x)$
  拉格朗日对偶函数：
  d ( λ , μ ) = m i n x ∈ X { f ( x ) + ∑ i = 1 m λ i g i ( x ) + ∑ j = 1 l μ j h j ( x ) } , 其中 X 为满足条件的 x 变量 ≤ m i n x ∈ S { f ( x ) + ∑ i = 1 m λ i g i ( x ) + ∑ j = 1 l μ j h j ( x ) } , 由于 g i ( x ) ≤ 0 , h j ( x ) = 0 , λ i ≥ 0 , 其中 S 为可行域 ≤ m i n x ∈ S { f ( x ) } d(\lambda,\mu) = min_{x\in X}\{ f(x) + \sum\limits_{i=1}^{m}\lambda_i g_i(x) + \sum\limits_{j=1}^{l}\mu_j h_j(x)\} ,其中X为满足条件的x变量\\ \le min_{x\in S}\{ f(x) + \sum\limits_{i=1}^{m}\lambda_i g_i(x) + \sum\limits_{j=1}^{l}\mu_j h_j(x) \},由于g_i(x) \le 0,h_j(x) = 0,\lambda_i \ge 0 ,其中S为可行域\\ \le min_{x\in S}\{f(x) \} d(λ,μ)=minx∈X​{f(x)+i=1∑m​λi​gi​(x)+j=1∑l​μj​hj​(x)},其中X为满足条件的x变量≤minx∈S​{f(x)+i=1∑m​λi​gi​(x)+j=1∑l​μj​hj​(x)},由于gi​(x)≤0,hj​(x)=0,λi​≥0,其中S为可行域≤minx∈S​{f(x)}
  因此：拉格朗日对偶函数$d(\lambda ,\mu )$是原问题最优解的函数值$p^{\ast }$的下界，即每个不同的$\lambda$与$\mu$确定的$d(\lambda ,\mu )$都是$p^{\ast }$的下界，但是我们希望下界越大越好，因为越大就更能接近真实的$p^{\ast }$。因此：
  拉格朗日对偶问题(D)转化为：
  $$\begin{gathered} ma{x}_{\lambda ,\mu }d(\lambda ,\mu ) \\ s.t.{\lambda }_i\ge 0,i=1,2,...,m \\ 也就是： \\ ma{x}_{\lambda \ge 0,\mu }mi{n}_{x\in S}L(x,\lambda ,\mu ) \end{gathered}$$
  我们可以观察到，对偶问题是关于$\lambda$和$\mu$的线性函数，因此对偶问题是一个凸优化问题，凸优化问题在最优化理论较为简单。
  弱对偶定理：对偶问题(D)的最优解$D^{\ast }$一定小于原问题最优解$P^{\ast }$，这点在刚刚的讨论得到了充分的证明，一定成立。
  强对偶定理：对偶问题(D)的最优解$D^{\ast }$在一定的条件下等于原问题最优解$P^{\ast }$，条件非常多样化且不是唯一的，也就是说这是个开放性的问题，在这里我给出一个最简单的条件，即：$f(x)$与$g_i(x)$为凸函数，$h_j(x)$为线性函数，X是凸集，$x^{\ast }$满足KKT条件，那么$D^{\ast }=P^{\ast }$。

  • 支持向量回归SVR
    在介绍完了相关的优化知识以后，我们开始正式学习支持向量回归SVR。
    

    在线性回归的理论中，每个样本点都要计算平方损失，但是SVR却是不一样的。SVR认为：落在$f(x)$的$ϵ$邻域空间中的样本点不需要计算损失，这些都是预测正确的，其余的落在$ϵ$邻域空间以外的样本才需要计算损失，因此：
    

$$\begin{gathered} mi{n}_{w,b,{\xi }_i,{\hat{\xi }}_i}\frac{1}{2}||w|{|}^2+C\sum _{i=1}^N({\xi }_i,{\hat{\xi }}_i) \\ s.t.f(x_i)-y_i\le ϵ+{\xi }_i \\ {y}_i-f(x_i)\le ϵ+{\hat{\xi }}_i \\ {\xi }_i,{\hat{\xi }}_i\le 0,i=1,2,...,N \end{gathered}$$

引入拉格朗日函数：
$$\begin{array}{l}L(w,b,\alpha ,\hat{\alpha },\xi ,\xi ,\mu ,\hat{\mu })\\ =\frac{1}{2}\Vert w{\Vert }^2+C{\sum }_{i=1}^N\left({\xi }_i+{\hat{\xi }}_i\right)-{\sum }_{i=1}^N{\xi }_i{\mu }_i-{\sum }_{i=1}^N{\hat{\xi }}_i{\hat{\mu }}_i\\ +{\sum }_{i=1}^N{\alpha }_i\left(f\left(x_i\right)-y_i-ϵ-{\xi }_i\right)+{\sum }_{i=1}^N{\hat{\alpha }}_i\left(y_i-f\left(x_i\right)-ϵ-{\hat{\xi }}_i\right)\end{array}$$
再令$L(w,b,\alpha ,\hat{\alpha },\xi ,\xi ,\mu ,\hat{\mu })$对$w,b,\xi ,\hat{\xi }$求偏导等于0，得： $w={\sum }_{i=1}^N\left({\hat{\alpha }}_i-{\alpha }_i\right)x_i$。
上述过程中需满足KKT条件，即要求：
$$\{\begin{array}{c}{\alpha }_i\left(f\left(x_i\right)-y_i-ϵ-{\xi }_i\right)=0\\ \widehat{{\alpha }_i}\left(y_i-f\left(x_i\right)-ϵ-{\hat{\xi }}_i\right)=0\\ {\alpha }_i{\hat{\alpha }}_i=0,{\xi }_i{\hat{\xi }}_i=0\\ \left(C-{\alpha }_i\right){\xi }_i=0,\left(C-{\hat{\alpha }}_i\right){\hat{\xi }}_i=0\end{array}$$
SVR的解形如:$f(x)={\sum }_{i=1}^N\left({\hat{\alpha }}_i-{\alpha }_i\right)x_i^{T}x+b$

sklearn中使用SVR实例：
sklearn.svm.SVR(*, kernel=‘rbf’, degree=3, gamma=‘scale’, coef0=0.0, tol=0.001, C=1.0, epsilon=0.1, shrinking=True, cache_size=200, verbose=False, max_iter=-1) Sklearn-SVR官网

• 参数：
  kernel：核函数，{‘linear’, ‘poly’, ‘rbf’, ‘sigmoid’, ‘precomputed’}, 默认=’rbf’。(后面会详细介绍)
  degree：多项式核函数的阶数。默认 = 3。
  C：正则化参数，默认=1.0。(后面会详细介绍)
  epsilon：SVR模型允许的不计算误差的邻域大小。默认0.1。

from sklearn.svm import SVR
from sklearn.preprocessing import StandardScaler     # 标准化数据
from sklearn.pipeline import make_pipeline   # 使用管道，把预处理和模型形成一个流程

reg_svr = make_pipeline(StandardScaler(), SVR(C=1.0, epsilon=0.2))
reg_svr.fit(X, y)
reg_svr.score(X,y)
#0.7024525421955277

2.5 优化基础模型

在刚刚的回归问题的基本算法中，我们使用数据集去估计模型的参数，如线性回归模型中的参数w，那么这个数据集我们称为训练数据集，简称训练集。我们在回归问题中使用训练集估计模型的参数的原则一般都是使得我们的损失函数在训练集达到最小值，其实在实际问题中我们是可以让损失函数在训练集最小化为0，如：在线性回归中，我加入非常多的高次项，使得我们模型在训练集的每一个数据点都恰好位于曲线上，那这时候模型在训练集的损失值也就是误差为0。

既然能做到这件事，是不是代表我们的建模完事大吉呢？换句话说我们的模型可以预测任意情况呢？答案是显然否定的。我们建立机器学习的目的并不是为了在已有的数据集，也就是训练集上效果表现非常优异，我们希望建立的机器学习模型在未知且情况复杂的测试数据上表现优异，我们称这样的未出现在训练集的未知数据集成为测试数据集，简称测试集。我们希望模型在测试集上表现优异！因为假如我们根据股票市场前六个月的数据拟合一个预测模型，我们的目的不是为了预测以前这六个月越准越好，而是预测明天乃至未来的股价变化。

2.5.1 训练均方误差与测试均方误差：

​ 在回归中，我们最常用的评价指标为均方误差，即：$MSE=\frac{1}{N}\sum _{i=1}^N(y_i-\hat{f}(x_i){)}^2$，其中$\hat{f}(x_i)$是样本$x_i$应用建立的模型$\hat{f}$预测的结果。如果我们所用的数据是训练集上的数据，那么这个误差为训练均方误差，如果我们使用测试集的数据计算的均方误差，我们称为测试均方误差。一般而言，我们并不关心模型在训练集上的训练均方误差，我们关心的是模型面对未知的样本集，即测试集上的测试误差，我们的目标是使得我们建立的模型在测试集上的测试误差最小。那我们如何选择一个测试误差最小的模型呢？这是个棘手的问题，因为在模型建立阶段，我们是不能得到测试数据的，比如：我们在模型未上线之前是不能拿到未知且真实的测试数据来验证我们的模型的。在这种情况下，为了简便起见，一些观点认为通过训练误差最小化来选择模型也是可行的。这种观点表面看上去是可行的，但是存在一个致命的缺点，那就是：一个模型的训练均方误差最小时，不能保证测试均方误差同时也很小。对于这种想法构造的模型，一般在训练误差达到最小时，测试均方误差一般很大！如图：

可以看到：当我们的模型的训练均方误差达到很小时，测试均方误差反而很大，但是我们寻找的最优的模型是测试均方误差达到最小时对应的模型，因此基于训练均方误差达到最小选择模型本质上是行不同的。正如上右图所示：模型在训练误差很小，但是测试均方误差很大时，我们称这种情况叫模型的过拟合。

2.5.2 偏差-方差的权衡：

从上图的测试均方误差曲线可以看到：测试均方误差曲线呈现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 )$为建模任务的难度，这个量在我们的任务确定后是无法改变的，也叫做不可约误差。那么模型的方差和偏差的平方和究竟是什么呢？所谓模型的方差就是：用不同的数据集去估计$f$时，估计函数的改变量。举个例子：我们想要建立一个线性回归模型，可以通过输入中国人身高去预测我们的体重。但是显然我们没有办法把全中国13亿人做一次人口普查，拿到13亿人的身高体重去建立模型。我们能做的就是从13亿中抽1000个样本进行建模，我们对这个抽样的过程重复100遍，就会得到100个1000人的样本集。我们使用线性回归模型估计参数就能得到100个线性回归模型。由于样本抽取具有随机性，我们得到的100个模型不可能参数完全一样，那么这100个模型之间的差异就叫做方差。显然，我们希望得到一个稳定的模型，也就是在不同的样本集估计的模型都不会相差太大，即要求f的方差越小越好。一般来说，模型的复杂度越高，f的方差就会越大。 如加入二次项的模型的方差比线性回归模型的方差要大。

​ 另一方面，模型的偏差是指：为了选择一个简单的模型去估计真实函数所带入的误差。假如真实的数据X与Y的关系是二次关系，但是我们选择了线性模型进行建模，那由于模型的复杂度引起的这种误差我们称为偏差，它的构成是复杂的。偏差度量了学习算法的期望预测与真实结果的偏离程度，即刻画了学习算法本身的拟合能力。偏差度量的是单个模型的学习能力，而方差度量的是同一个模型在不同数据集上的稳定性。“偏差-方差分解”说明：泛化性能是由学习算法的能力、数据的充分性以及学习任务本身的难度所共同决定的。给定学习任务，为了取得好的泛化性能，则需使偏差较小，即能够充分拟合数据，并且使方差较小，即使得数据扰动产生的影响小。

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

2.5.3 特征提取：

​ 在前面的讨论中，我们已经明确一个目标，就是：我们要选择一个测试误差达到最小的模型。但是实际上我们很难对实际的测试误差做精确的计算，因此我们要对测试误差进行估计，估计的方式有两种：训练误差修正与交叉验证。

$MSE=\frac{1}{N}\sum _{i=1}^N(y_i-\hat{f}(x_i){)}^2$
$$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)$的偏差平方和误差项$ϵ$的方差

训练误差修正（间接估计）：

前面的讨论我们已经知道，模型越复杂，训练误差越小，测试误差先减后增。因此，我们先构造一个特征较多的模型使其过拟合，此时训练误差很小而测试误差很大，那这时我们加入关于特征个数的惩罚。因此，当我们的训练误差随着特征个数的增加而减少时，惩罚项因为特征数量的增加而增大，抑制了训练误差随着特征个数的增加而无休止地减小。有几种间接估计测试误差的方法：

1. $C_p=\frac{1}{N}(RSS+2d{\hat{\sigma }}^2)$

其中d为模型特征个数，$RSS=\sum _{i=1}^N(y_i-\hat{f}(x_i){)}^2$$RSS$ 也可以理解为训练均方误差MSE。${\hat{\sigma }}^2$为模型预测误差的方差的估计值，即残差的方差。$2d{\hat{\sigma }}^2$ 为惩罚项。当模型特征个数增加时，训练误差减小，即RSS减小，但d增大，是整个$C_p$靠近测试误差

2. AIC赤池信息量准则：$AIC=\frac{1}{d{\hat{\sigma }}^2}(RSS+2d{\hat{\sigma }}^2)$
3. BIC贝叶斯信息量准则：$BIC=\frac{1}{n}(RSS+log(n)d{\hat{\sigma }}^2)$

$C_p$和AIC、BIC都是作为训练误差修正系列的具体方法用来间接估计测试误差，BIC是惩罚最终的，最终会得到特征最少的模型。间接估计简单的说就是对训练误差进行修正来靠近真实的测试误差。

交叉验证（直接估计）：

前面讨论的对训练误差修正得到测试误差的估计是间接方法，这种方法的桥梁是训练误差，而交叉验证则是对测试误差的直接估计。交叉验证比训练误差修正的优势在于：能够给出测试误差的一个直接估计。在这里只介绍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(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​}中测试误差最小的模型作为最优模型。

2.5.4 压缩估计(正则化Regularization)：

​ 除了刚刚讨论的直接对特征自身进行选择以外，我们还可以对回归的系数进行约束或者加罚的技巧对p个特征的模型进行拟合，显著降低模型方差，这样也会提高模型的拟合效果。具体来说，就是将回归系数往零的方向压缩（类似特征提取），这也就是为什么叫压缩估计的原因了。（一些特征的系数为0，所以该特征近似于被忽略）

2.5.4.1 岭回归(L2正则化)：

在线性回归中，我们的损失函数为$J(w)=\sum _{i=1}^N(y_i-w_0-\sum _{j=1}^p{w}_j{x}_{ij}{)}^2$，我们在线性回归的损失函数的基础上添加对系数的约束或者惩罚（L2范式），即：
$$\begin{gathered} \hat{w}=(X^{T}X{)}^{-1}{X}^{T}Y,线性回归中参数的估计，当样本数小于特征维度时不能求逆 \\ 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，I为单位矩阵 \end{gathered}$$

​ 调节参数$\lambda$的大小是影响压缩估计的关键，$\lambda$越大，惩罚的力度越大，系数则越趋近于0，反之，选择合适的$\lambda$对模型精度来说十分重要。岭回归通过牺牲线性回归的无偏性降低方差（牺牲偏差，降低方差，且方差减小的量大于偏差增大的量），有可能使得模型整体的测试误差较小，提高模型的泛化能力。

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

​ 最优化参数用到LARS最小角回归。岭回归的一个很显著的特点是：将模型的系数往零的方向压缩，但是岭回归的系数只能呢个趋于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约束(L1为绝对值约束，所以是矩形；L2为平方项约束，所以是圆形)，Lsaao回归和岭回归都是在约束下的回归，因此最优的参数为椭圆形曲线（红色椭圆为RSS-训练误差的等高线）与菱形和圆形区域相切的点。但是Lasso回归的约束在每个坐标轴上都有拐角，因此当RSS曲线与坐标轴相交时恰好回归系数中的某一个为0，这样就实现了特征提取。反观岭回归的约束是一个圆域，没有尖点，因此与RSS曲线相交的地方一般不会出现在坐标轴上，因此无法让某个特征的系数为0，因此无法做到特征提取。

2.5.5 降维：

​ 到目前为止，我们所讨论的方法对方差的控制有两种方式：一种是使用原始变量的子集，另一种是将变量系数压缩至零。但是这些方法都是基于原始特征$x_1,...,x_p$得到的，现在我们探讨一类新的方法：将原始的特征空间投影到一个低维的空间实现变量的数量变少，如：将二维的平面投影至一维空间。机器学习领域中所谓的降维就是指采用某种映射方法，将原高维空间中的数据点映射到低维度的空间中。降维的本质是学习一个映射函数 f : x->y，其中x是原始数据点的表达，目前最多使用向量表达形式。 y是数据点映射后的低维向量表达，通常y的维度小于x的维度（当然提高维度也是可以的）。f可能是显式的或隐式的、线性的或非线性的。目前大部分降维算法处理向量表达的数据，也有一些降维算法处理高阶张量表达的数据。之所以使用降维后的数据表示是因为在原始的高维空间中，包含有冗余信息以及噪音信息，在实际应用例如图像识别中造成了误差，降低了准确率；而通过降维,我们希望减少 冗余信息 所造成的误差,提高识别（或其他应用）的精度。又或者希望通过降维算法来寻找数据内部的本质结构特征。在很多算法中，降维算法成为了数据预处理的一部分，如PCA。事实上，有一些算法如果没有降维预处理，其实是很难得到很好的效果的。 (摘自：rosenor1博客)

2.5.5.1 主成分分析(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个特征向量即可。

特征提取的实例：

向前逐步回归实例分享：

案例来源：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.api as sm #最小二乘
from statsmodels.formula.api import ols #加载ols模型
forward_select(data=boston_data,target="Price")

​

​

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

岭回归实例分享：

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.739957023371629

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.7140164719858566

(5) 对模型超参数进行调优(调参)：
在刚刚的讨论中，我们似乎对模型的优化都是对模型算法本身的改进，比如：岭回归对线性回归的优化在于在线性回归的损失函数中加入L2正则化项从而牺牲无偏性降低方差。但是，大家是否想过这样的问题：在L2正则化中参数$\lambda$应该选择多少？是0.01、0.1、还是1？到目前为止，我们只能凭经验或者瞎猜，能不能找到一种方法找到最优的参数$\lambda$？事实上，找到最佳参数的问题本质上属于最优化的内容，因为从一个参数集合中找到最佳的值本身就是最优化的任务之一，我们脑海中浮现出来的算法无非就是：梯度下降法、牛顿法等无约束优化算法或者约束优化算法，但是在具体验证这个想法是否可行之前，我们必须先认识两个最本质概念的区别。

• 参数与超参数：
  我们很自然的问题就是岭回归中的参数$\lambda$和参数w之间有什么不一样？事实上，参数w是我们通过设定某一个具体的$\lambda$后使用类似于最小二乘法、梯度下降法等方式优化出来的，我们总是设定了$\lambda$是多少后才优化出来的参数w。因此，类似于参数w一样，使用最小二乘法或者梯度下降法等最优化算法优化出来的数我们称为参数，类似于$\lambda$一样，我们无法使用最小二乘法或者梯度下降法等最优化算法优化出来的数我们称为超参数。
  模型参数是模型内部的配置变量，其值可以根据数据进行估计。
  • 进行预测时需要参数。
  • 它参数定义了可使用的模型。
  • 参数是从数据估计或获悉的。
  • 参数通常不由编程者手动设置。
  • 参数通常被保存为学习模型的一部分。
  • 参数是机器学习算法的关键，它们通常由过去的训练数据中总结得出 。
    模型超参数是模型外部的配置，其值无法从数据中估计。
  • 超参数通常用于帮助估计模型参数。
  • 超参数通常由人工指定。
  • 超参数通常可以使用启发式设置。
  • 超参数经常被调整为给定的预测建模问题。
    我们前面(4)部分的优化都是基于模型本身的具体形式的优化，那本次(5)调整的内容是超参数，也就是取不同的超参数的值对于模型的性能有不同的影响。
• 网格搜索GridSearchCV()：
  网格搜索：https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html?highlight=gridsearchcv#sklearn.model_selection.GridSearchCV
  网格搜索结合管道：https://scikit-learn.org/stable/auto_examples/compose/plot_compare_reduction.html?highlight=gridsearchcv
  网格搜索的思想非常简单，比如你有2个超参数需要去选择，那你就把所有的超参数选择列出来分别做排列组合。举个例子：$\lambda =0.01,0.1,1.0$和$\alpha =0.01,0.1,1.0$,你可以做一个排列组合，即：{[0.01,0.01],[0.01,0.1],[0.01,1],[0.1,0.01],[0.1,0.1],[0.1,1.0],[1,0.01],[1,0.1],[1,1]} ，然后针对每组超参数分别建立一个模型，然后选择测试误差最小的那组超参数。换句话说，我们需要从超参数空间中寻找最优的超参数，很像一个网格中找到一个最优的节点，因此叫网格搜索。
• 随机搜索 RandomizedSearchCV() ：
  https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.RandomizedSearchCV.html?highlight=randomizedsearchcv#sklearn.model_selection.RandomizedSearchCV
  网格搜索相当于暴力地从参数空间中每个都尝试一遍，然后选择最优的那组参数，这样的方法显然是不够高效的，因为随着参数类别个数的增加，需要尝试的次数呈指数级增长。有没有一种更加高效的调优方式呢？那就是使用随机搜索的方式，这种方式不仅仅高效，而且实验证明，随机搜索法结果比稀疏化网格法稍好(有时候也会极差，需要权衡)。参数的随机搜索中的每个参数都是从可能的参数值的分布中采样的。与网格搜索相比，这有两个主要优点：
  • 可以独立于参数数量和可能的值来选择计算成本。
  • 添加不影响性能的参数不会降低效率。

SVR调优

# 我们先来对未调参的SVR进行评价： 
import numpy as np
from sklearn.svm import SVR     # 引入SVR类
from sklearn.pipeline import make_pipeline   # 引入管道简化学习流程
from sklearn.preprocessing import StandardScaler # 由于SVR基于距离计算，引入对数据进行标准化的类
from sklearn.model_selection import GridSearchCV  # 引入网格搜索调优
from sklearn.model_selection import cross_val_score # 引入K折交叉验证
from sklearn import datasets

boston = datasets.load_boston()     # 返回一个类似于字典的类
X = boston.data
y = boston.target
features = boston.feature_names
pipe_SVR = make_pipeline(StandardScaler(),SVR())
score1 = cross_val_score(estimator=pipe_SVR,X = X,y = y,scoring = 'r2',cv = 10)       # 10折交叉验证
print("CV accuracy: %.3f +/- %.3f" % ((np.mean(score1)),np.std(score1)))
#CV accuracy: 0.187 +/- 0.649

网格搜索来对SVR调参

# 下面我们使用网格搜索来对SVR调参：
from sklearn.pipeline import Pipeline
pipe_svr = Pipeline([("StandardScaler",StandardScaler()),("svr",SVR())])
param_range = [0.0001,0.001,0.01,0.1,1.0,10.0,100.0,1000.0]
param_grid = [{"svr__C":param_range,"svr__kernel":["linear"]},{"svr__C":param_range,"svr__gamma":param_range,"svr__kernel":["rbf"]}]# 注意__是指两个下划线，一个下划线会报错的
gs = GridSearchCV(estimator=pipe_svr,param_grid = param_grid,scoring = 'r2',cv = 10)       # 10折交叉验证
gs = gs.fit(X,y)
print("网格搜索最优得分：",gs.best_score_)
print("网格搜索最优参数组合：\n",gs.best_params_)
#网格搜索最优得分： 0.6081303070817233
#网格搜索最优参数组合：{'svr__C': 1000.0, 'svr__gamma': 0.001, 'svr__kernel': 'rbf'}

随机搜索来对SVR调参：

# 下面我们使用随机搜索来对SVR调参：
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import uniform  # 引入均匀分布设置参数
pipe_svr = Pipeline([("StandardScaler",StandardScaler()),
                                                         ("svr",SVR())])
distributions = dict(svr__C=uniform(loc=1.0, scale=4),    # 构建连续参数的分布
                     svr__kernel=["linear","rbf"],                                   # 离散参数的集合
                    svr__gamma=uniform(loc=0, scale=4))

rs = RandomizedSearchCV(estimator=pipe_svr,
                                                     param_distributions = distributions,
                                                     scoring = 'r2',
                                                      cv = 10)       # 10折交叉验证
rs = rs.fit(X,y)
print("随机搜索最优得分：",rs.best_score_)
print("随机搜索最优参数组合：\n",rs.best_params_)
#随机搜索最优得分： 0.30021249798866756
#随机搜索最优参数组合：{'svr__C': 1.4195029566223933, 'svr__gamma': 1.8683733769303625, 'svr__kernel': 'linear'}