【学习笔记】集成学习（三）：偏差和方差理论

Datawhale组队学习第27期：集成学习
本次学习的指导老师萌弟的教学视频
本贴为学习记录帖，有任何问题欢迎随时交流~
部分内容可能还不完整，后期随着知识积累逐步完善。
开始时间：2021年7月13日
最新更新：2021年7月19日（Task3偏差和方差理论）

一、训练误差与测试误差

偏差和方差理论

• 偏差-方差的权衡
  $$E(y_0-\hat{f}(x_0){)}^2=Var(\hat{f}(x_0))+[Bias(\hat{f}(x_0)){]}^2+Var(ϵ)$$

• 分解为$\hat{f}(x_0)$的方差、其偏差平方和和误差项。

• 测试均方误差的期望不会低于误差的方差

  • $Var(ϵ)$是建模任务的难度，不可约误差
    • 抽出样本集之间的误差是模型的方差（Var），越小越好

• 模型复杂度与模型的方差正相关，与模型的偏差负相关。

• 模型复杂度引起的误差称为偏差（Bias），提高泛化能力则需要使得偏差变小

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

• 泛化性能是由学习算法的能力（偏差）、数据的充分性（方差）和学习任务本身难度（不可改变）共同决定的。

二、模型简化

特征提取

• 测试误差的估计

  • 直接估计（交叉验证）
  • 间接估计（训练误差修正）
    • 模型越复杂，训练误差越小，测试误差会先增后减
    • 先过拟合，再加入特征个数的惩罚
    • 如，$C_p$、$AIC$、$BIC$

• 从p个特征中选择m个特征，使得模型的测试误差的估计最小

  • 最优子集选择
    • 确定基模型，记为$M_0$
    • 逐次增加一个变量，每次记录增加一个变量时测试误差最小的那个模型，记为$M_k$
    • 直到增加了p个，最后共有$p+1$个模型，再进行比较，选择测试误差最小的那个。
    • 最终总共需要$2^p$次，计算量和内存消耗大。
  • 向前逐步选择
    • 确定基模型，记为$M_0$
    • 逐次增加一个变量，而每次增加的变量是在上一个测试误差最小的模型上增加，与最优子集不同。
    • 总共需要$1+2+...+p$次

压缩估计（正则化）

• 训练中加入约束（惩罚项），回归系数往零的方向压缩

• 岭回归

  • L2正则化
  • 使得w计算中可逆求解
  • 牺牲线性回归中的无偏性来降低方差，有可能使得模型整体的测试误差较小，提高模型的泛化能力
  • 无法做特征选择？

• Lasso回归

  • L1正则化

  • 不可以求导

  • LARS最小角回归

• 在二维的情况下，两者区别如下：

降维

• 高维度空间映射到低维度的空间

• 一种是降噪，提高识别精度；一种是寻找数据内部的本质结构特征

• PCA（主成分）

  • 中心化
  • 计算每个点在某个方向上的投影
  • 计算投影方差
  • 最大化投影方差，求出对应的方向$u_1$

  $$\begin{gathered} S^2{u}_1=\lambda u_1 \\ \lambda 是特征值，u_1是特征向量 \end{gathered}$$

  • 若特征有m个，取p个特征值，从而达到降维的效果

三、调参

• 如，L2正则化中有超参数$\lambda$

• 超参数预先需要设定的值，而非优化问题的参数值

• 超参数是人工指定、帮助估计模型参数，启发式设置等

• 常用方法

  • 网格搜索：$GridSearchCV()$
  • 随机搜索：$RandomizedSearchCV()$
  • 贝叶斯搜索

四、Statsmodels包的使用

• statsmodels的OLS使用方法：

  • 包中主要有两个模块可以实现线性回归：
    • statsmodels.api.OLS
    • statsmodels.formula.api.ols

• 对于statsmodels.api.OLS

  • 该模块可以直接传入自变量X和因变量Y，然后求值计算。
  • 该用法结果不带截距项，加入截距项的话需要先对数据集使用statsmodels.api.add_constant()增加额外一列。

• statsmodels.formula.api.ols
  该模块是根据变量名（formula）进行传值计算的，传入的数据集有要求。
  

  该用法是带有截距项的，如果需要取消截距项，在formula中加入-1，从而排除截距项。

  # 两种等价写法（前者自动加入截距项，后者需要手动增加）
  formula = '{}~{} + {}'.format(d_set.columns[-1], d_set.columns[0], d_set.columns[1])
  res = ols(formula, data=d_set).fit()
  res.summary()
  
  # sm.add_constant(d_x.values[:, :2])这是一种加入截距项的方法
  res2 = sm.OLS(d_y.values, d_x.values[:, :2]).fit()
  res2.summary()
  

五、关于测试误差的间接估计

Log Likelihood

• 这里首先回顾线性回归的基本假设。

  • 线性假定：$y=w^{T}x+ϵ$，其中w的元素均为常数

  • 严格外生性假定：$E(ϵ|x)=0$，即扰动项独立于观测数据

  • 数值矩阵 $x$ 是列满秩，不存在严格的多重共线性

  • 球形扰动项假定：$Var(ϵ|x)={\sigma }^2$，即扰动项条件同方差和无自相关

  • 正态性假定：$ϵ|x\sim N(0,{\sigma }^2)$

• 然后由线性回归模型 $y=w^{T}x+ϵ$ 可知，$y|x,w\sim N(w^{T}x,{\sigma }^2)$，根据极大似然估计和$\hat{y}=w^{T}x$可知：
  ln ⁡ L ( y ∣ x ; w ) = ln ⁡ P ( y ∣ x ; w ) = ln ⁡ ∏ i = 1 n P ( y i ∣ x i ; w ) = ∑ i = 1 n ln ⁡ ( 1 2 π σ 2 e x p { − ( y i − w T x i ) 2 2 σ 2 } ) = ∑ i = 1 n ln ⁡ ( 1 2 π σ 2 − ( y i − w T x i ) 2 2 σ 2 ) = − n 2 ln ⁡ ( 2 π σ 2 ) − ∑ i = 1 n ( y i − y ^ i ) 2 2 σ 2 = − n 2 ln ⁡ ( 2 π ) − n 2 ln ⁡ ( σ 2 ) − 1 2 σ 2 ∑ i = 1 n ( y i − y ^ i ) 2 \begin{aligned} \ln L(y|x; w) &= \ln P(y|x;w) \\ &= \ln \prod_{i=1}^{n}P(y_i|x_i;w) \\ &= \sum\limits_{i=1}^{n}\ln(\frac{1}{\sqrt{2\pi}\sigma^2} exp\{- \frac{(y_i - w^Tx_i)^2}{2\sigma^2} \}) \\ &= \sum\limits_{i=1}^{n}\ln(\frac{1}{\sqrt{2\pi}\sigma^2} - \frac{(y_i - w^Tx_i)^2}{2\sigma^2}) \\ &= - \frac{n}{2}\ln(2\pi \sigma^2) - \sum\limits_{i=1}^{n} \frac{(y_i-\hat y_i)^2}{2\sigma^2} \\ &= - \frac{n}{2}\ln(2\pi) - \frac{n}{2}\ln(\sigma^2) - \frac{1}{2\sigma^2} \sum\limits_{i=1}^{n}(y_i-\hat y_i)^2 \end{aligned} lnL(y∣x;w)​=lnP(y∣x;w)=lni=1∏n​P(yi​∣xi​;w)=i=1∑n​ln(2π ​σ21​exp{−2σ2(yi​−wTxi​)2​})=i=1∑n​ln(2π ​σ21​−2σ2(yi​−wTxi​)2​)=−2n​ln(2πσ2)−i=1∑n​2σ2(yi​−y^​i​)2​=−2n​ln(2π)−2n​ln(σ2)−2σ21​i=1∑n​(yi​−y^​i​)2​

• 实际上，$w$可以通过转化为最优化问题求值。除此之外，${\sigma }^2$也是需要估计的参数，一般可以通过代入参数$w$进行估计。

  $$\frac{\partial lnL}{\partial {\sigma }^2}=-\frac{n}{2{\sigma }^2}+\frac{1}{2({\sigma }^2{)}^2}\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2=0$$

• 解得${\sigma }^2$的极大似然估计为：

  $${\hat{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$$

• 将${\hat{\sigma }}^2$和$\hat{y}$代入对数似然函数中，可得：

  $$\begin{array}{rl}\ln L(y|x;w)& =-\frac{n}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(y_i-\hat{y}{)}^2\\ & =-\frac{n}{2}\ln (2\pi {\sigma }^2)-\frac{n}{2{\sigma }^2}\cdot \frac{1}{n}\sum _{i=1}^n(y_i-\hat{y}{)}^2\\ & =-\frac{n}{2}\ln (2\pi {\sigma }^2)-\frac{n}{2{\sigma }^2}\cdot {\sigma }^2\\ & =-\frac{n}{2}\ln (2\pi {\sigma }^2)-\frac{n}{2}\\ & =-\frac{n}{2}\ln (2\pi )-\frac{n}{2}\ln ({\sigma }^2)-\frac{n}{2}\end{array}$$

• 显然，记 $n{\sigma }^2$ 为$SSE$ (误差平方和，也有的称为$RSS$)，则有：
  $$\ln L(\mu ,{\sigma }^2|X)=-\frac{n}{2}\ln (2\pi )-\frac{n}{2}\ln (\frac{SSR}{n})-\frac{n}{2}$$

• 事实上，由极大似然估计的思想来看，模型比较时也可以使用Log Likelihood（值越大，模型越好），统计软件一般会给出该值。

AIC

• 赤池信息准则$AIC$的计算方法：
  $$\begin{array}{rl}AIC& =2k-2ln(L),其中k为参数的个数\\ \ln (L)& =-\frac{n}{2}ln(2\pi )-\frac{n}{2}ln(\frac{SSE}{n})-\frac{n}{2}=-\frac{n}{2}[ln(2\pi )+ln(\frac{SSE}{n})+1]\\ SSE& =\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2\end{array}$$

六、向前逐步回归的实现

注意点

• 开源项目中选择$AIC$作为模型优良性的判断标准，这里参照上面的知识点，$AIC$是间接估计，还可以采用交叉验证法对模型的测试误差进行直接估计。

• 开源项目中给出的代码是使用statsmodels.formula.api.ols，本文给出的是statsmodels.api.OLS的方法

原理

• 初始化模型$M_0$，计算此时的测试误差，可以考虑设置为无穷大，即：float('inf')​
• 初次先加入单个变量拟合模型，找到使得模型的$AIC$值最小的变量，作为模型$M_1$
• 在$M_1$模型的基础上，再加入单个新变量，找到使得模型的$AIC$值最小的变量，作为模型$M_2$。
• 按上述步骤遍历所有遍历，直到加入所有变量（假定总共p个特征变量），最终得到 { M 1 , M 2 , . . . , M p } \{M_1, M_2, ..., M_p\} {M1​,M2​,...,Mp​}共p个模型。
• 最后找到这p个模型中$AIC$值达到最小的模型。

算法设计

• 输入：特征变量X，标签Y，特征变量名feature_name

• 初始化：总模型得分md_scores、变量遍历feature、变量选择choose

• 计算$M_0$的模型得分（只有截距项，截距项为 $y$ 均值）

• 设置外层循环，终止条件为：feature为空

  • 初始化：单轮模型得分scores
  • 设置内层循环，遍历md_feature
    • 每次选定一个变量加入模型中，模型得分和变量名加入到scores中
  • 寻找模型得分的最大值加入到md_scores中，对应的变量名加入到choose中。
  • feature弹出选定的变量。

• 遍历md_scores，找到最大值所在的下标，输出choose到对应下标的元素列表。

代码实现

import numpy as np
import pandas as pd
from sklearn.datasets import load_boston
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 载入数据集
data = load_boston()

# 设定随机种子（后面用于生成随机参数）
np.random.seed(2021)
np.set_printoptions(suppress=True)

# 选取变量
x = data['data']
y = data['target']
d_x = pd.DataFrame(x, columns=data['feature_names'])
d_y = pd.DataFrame(y, columns=['PRICE'])
d_set = pd.concat([d_x, d_y], axis=1)

# 定义向前逐步回归（加入截距项）
def forward_reg2(x_, y_, cons=False):
    md_scores = list()
    choose = list()
    train_x = x_
    if cons:
        train_x = sm.add_constant(x_, prepend=True)
        choose.append('const')
    feature = list(x_.keys())
    md_scores.append(sm.OLS(y_, np.ones(len(y))).fit().aic)     # 模型不含特征变量时
    loop = 0
    while feature:
        best_score = float('inf')   # 初始化为无穷大
        choose_idx = 0   # 初始化单轮选择的变量名
        for i in range(len(feature)):
            score = sm.OLS(y_, train_x[choose+[feature[i]]]).fit().aic
            if score < best_score:
                best_score = score
                choose_idx = i
        choose.append(feature.pop(choose_idx))
        md_scores.append(best_score)
        loop += 1
        print('第{}轮，加入变量{}，AIC为：{}'.format(loop, choose[-1], best_score))
    best_score = float('inf')
    best_idx = 0
    for idx, md_score in enumerate(md_scores):
        if md_score < best_score:
            best_idx = idx
            best_score = md_score
    print('最终选取特征变量为：{}'.format(choose[1:best_idx +1]))
    return choose[1:best_idx + 1]


ols_feature = forward_reg2(d_x, d_y, True)

• 最终结果如下：
  

七、多个模型比较

• 上述提到了模型的简化主要有：

  • 特征提取，通过对测试误差进行估计，筛选变量。如，最优子集选择
  • 压缩估计，对损失函数加入惩罚项。如，岭回归和Lasso回归
  • 降维，高维度空间映射到低维度的空间。如，PCA降维

• 为了更好得比较三种简化方法，实验加入了以下工具：

  • 利用train_test_split进行分割训练集和测试集，测试集用于评估模型的精度。（注：这里暂不考虑调参的因素，并未考虑验证集）
  • 利用pipeline对向前逐步回归进行重新定义，定义方法参考该文章
  • 模型的评价指标统一采用sklearn.metrics.r2_score

• 完整代码如下所示：

  import numpy as np
  import pandas as pd
  from sklearn.datasets import load_boston
  import os
  from sklearn.linear_model import LinearRegression, Lasso, Ridge
  from sklearn.decomposition import PCA
  from sklearn.preprocessing import StandardScaler
  from sklearn.pipeline import make_pipeline
  from sklearn.model_selection import train_test_split
  import statsmodels.api as sm
  from sklearn.base import BaseEstimator, RegressorMixin
  from sklearn.metrics import r2_score
  
  
  # 定义向前逐步回归
  class ForwardPipeline(BaseEstimator, RegressorMixin):
      def __init__(self):
          """
          初始化，参考BaseEstimator的定义方式
          """
          # self.feature = None
          self.md = None
          self.transform = None
          self.best_feature = None
          self.mean = 0
          self.std = 0
  
      def fit(self, x_, y_):
          """
          拟合最优模型，注释的文字是可以进行改进的思路
          """
          # 记录变量标签名
          # self.feature = list(x_.keys())
  
          # 数据标准化
          # self.mean = np.mean(x_, axis=0)
          # self.std = np.std(x_, axis=0)
          # std_x = (x_ - self.mean) / self.std
          # std_x = pd.DataFrame(std_x, columns=self.feature)
          std_x = x_
  
          # 求最优参数
          self.best_feature = self.forward_reg(std_x, y_)
          self.md = sm.OLS(y_, std_x[self.best_feature]).fit()
          return self
  
      def predict(self, x_):
          # std_x = pd.DataFrame(self.transform.transform(x_), columns=self.feature)
          # std_x = (x_ - self.mean) / self.std
          std_x = x_
          y_pred = self.md.predict(std_x[self.best_feature])
          return y_pred
  
      def score(self, x_, y_, sample_weight=None):
          """
          计算R2，参考RegressorMixin的定义方式，对RegressorMixin的方法进行重写
          """
          y_pred = self.predict(x_)
          return r2_score(y_, y_pred, sample_weight=sample_weight)
  
      @staticmethod
      def forward_reg(x_, y_, cons=True):
          md_scores = list()
          choose = list()
          train_x = x_
          if cons:
              train_x = sm.add_constant(x_, prepend=True)
              choose.append('const')
          feature = list(x_.keys())
          md_scores.append(sm.OLS(y_, np.ones(len(y_))).fit().aic)  # 模型不含特征变量时
          loop = 0
          while feature:
              best_score = float('inf')  # 初始化为无穷大
              choose_idx = 0  # 初始化单轮选择的变量名
              for i in range(len(feature)):
                  score = sm.OLS(y_, train_x[choose + [feature[i]]]).fit().aic
                  if score < best_score:
                      best_score = score
                      choose_idx = i
              choose.append(feature.pop(choose_idx))
              md_scores.append(best_score)
              loop += 1
              # print('第{}轮，加入变量{}，AIC为：{}'.format(loop, choose[-1], best_score))
          best_score = float('inf')
          best_idx = 0
          for idx, md_score in enumerate(md_scores):
              if md_score < best_score:
                  best_idx = idx
                  best_score = md_score
          # print('最终选取特征变量为：{}'.format(choose[1:best_idx + 1]))
          return choose[1:best_idx + 1]
  
  
  if __name__ == '__main__':		
  	# 设定项目路径
  	os.chdir('D:/NewProject/datawhale/ensemble/task3')
  	np.set_printoptions(suppress=True)
  	np.random.seed(2021)
  	
  	# 加载数据集
  	data = load_boston()
  	print(data.keys())  # 查看数据集的字段
  	
  	# 选定变量
  	x = data['data']  # 返回的是ndarray类型
  	y = data['target']
  	
  	# 转化成dataframe结构
  	d_x = pd.DataFrame(x, columns=data['feature_names'])
  	d_y = pd.DataFrame(y, columns=['PRICE'])
  	data_set = pd.concat([d_x, d_y], axis=1)
  	
  	# 分割数据集
  	x_train, x_test, y_train, y_test = train_test_split(d_x, d_y, test_size=0.2)
  	
  	# 普通ols线性回归模型
  	md1 = make_pipeline(LinearRegression())
  	md1.fit(x_train, y_train)
  	print('普通ols回归的R2：', md1.score(x_test, y_test))
  	
  	# 岭回归，超参数alpha（又称为岭迹）
  	md2 = make_pipeline(StandardScaler(), Ridge(alpha=0.5))
  	md2.fit(x_train, y_train)
  	print('岭回归的R2：', md2.score(x_test, y_test))
  	
  	# Lasso回归，超参数alpha（正则化系数）
  	md3 = make_pipeline(StandardScaler(), Lasso(alpha=0.5))
  	md3.fit(x_train, y_train)
  	print('Lasso回归的R2：', md3.score(x_test, y_test))
  	
  	# PCA降维后回归，这里降维后选取的主成分个数也算是超参数
  	md4 = make_pipeline(StandardScaler(), PCA(n_components=0.8), LinearRegression())
  	md4.fit(x_train, y_train)
  	print('PCA降维后回归的R2：', md4.score(x_test, y_test))
  
  	# 向前逐步回归，这里并未带有超参数
  	md5 = make_pipeline(ForwardPipeline())
  	md5.fit(x_train, y_train)
  	print('向前逐步回归的R2：', md5.score(x_test, y_test))
  

• 显示结果如下：

  • 普通OLS效果低于岭回归和Lasso回归，而要高于PCA和向前逐步回归。
  • 一定程度上来说，PCA和向前逐步回归容易丢失原始数据的信息，从而降低了模型的精度。
  • 加入了正则化项后，模型精度有所提升，在后续进行调参时主要考虑前三个模型的比较。
    

八、线性回归模型调参

模型调整思路

• 上述模型的参数中如Lasso回归中的正则化系数$\alpha$是超参数，最开始是由我们自行设定的而非优化问题中求解的对象，而不同的超参数会影响模型的效果。
• 上述模型用所有的训练集和设定的参数直接得出模型，实际上是不符合研究规范的，具体可以看“如何正确使用机器学习中的训练集、验证集和测试集？”这篇文章了解训练集、验证集和测试集三者在机器学习中的作用。
• 模型进行超参数优化时需要充分使用验证集，通常采用交叉验证去得到模型的最优参数。

超参数优化

• 现只考虑岭回归和Lasso回归两个模型。
  • 岭回归模型中，超参数主要为：$alpha$，通常称为岭迹。
  • Lasso模型中，超参数主要为：$alpha$，称为正则化系数。
• 超参数优化实现包
  • 网格搜索GridSearchCV
  • 针对上述模型特有的超参数优化包：LassoCV和RidgeCV
• 简单实现代码：

# 自定义分割训练集和验证集的方法
kf = ShuffleSplit(n_splits=10, test_size=0.2, random_state=2021)	# 随机划分

# 对Lasso回归进行调参
lasso = Lasso()
lasso_param = {'alpha': [i for i in np.arange(0.1, 1, 0.01)]}
lasso_gcv = GridSearchCV(lasso, lasso_param, n_jobs=-1, cv=kf, scoring='r2', verbose=1)
lasso_gcv_res = lasso_gcv.fit(x_train, y_train)
print('Lasso回归最优参数：', lasso_gcv_res.best_params_)
print('Lasso回归得分：', lasso_gcv_res.score(x_test, y_test))

# 对岭回归进行调参
ridge = Ridge()
ridge_param = {'alpha': [i for i in np.arange(0.1, 1, 0.05)]}
ridge_gcv = GridSearchCV(ridge, ridge_param, n_jobs=-1, cv=kf, scoring='r2', verbose=1)
ridge_gcv_res = ridge_gcv.fit(x_train, y_train)
print('岭回归回归最优参数：', ridge_gcv_res.best_params_)
print('岭回归得分：', ridge_gcv_res.score(x_test, y_test))

# 基于两个模型的专用包
ls_ = LassoCV(eps=0.0001, n_alphas=300, cv=kf).fit(x_train, y_train)
ridge_ = RidgeCV(alphas=np.logspace(-4, 4, 100), cv=kf).fit(x_train, y_train)

• 训练的结果如下：

调参可视化

• 该部分是改编自sklearn的官方文档，用于展现参数与精度关系。
• 基于LassoCV，官方文档的命令（使用了matplotlib.pyplot）

  epsilon = 1e-4  # 用于避免semilogx中x为0
  plt.figure()
  plt.semilogx(ls_.alphas_ + epsilon, ls_.mse_path_, ':')	# 折线图，指数型（CV）
  plt.plot(ls_.alphas_ + epsilon, ls_.mse_path_.mean(axis=-1), 'k',
           label='Average across the folds', linewidth=2)	# 普通线图（mean）
  plt.axvline(ls_.alpha_ + epsilon, linestyle='--', color='k',
              label='alpha: CV estimate')	# 竖线图，最优参数所在位置
  plt.legend()
  plt.xlabel(r'$\alpha$')
  plt.ylabel('Mean square error')
  plt.title('Mean square error on each fold: coordinate descent ')
  plt.axis('tight')
  plt.show()
  

• 对上述代码重写，自定义可视化函数（这里的精度是r2）

  # 定义可视化函数
  def plot_cv(estimator, x_, y_, params_, cv, n_jobs=-1, scoring='r2', n_sp=10):
      md = estimator()	# 模型初始化
      md_gcv = GridSearchCV(md, params_, cv=cv, n_jobs=n_jobs, scoring=scoring, verbose=0)	# 调参初始化
      md_gcv_res = md_gcv.fit(x_, y_)	 # 寻找最优超参数
      alphas_ = params_['alpha']	# 参数空间
      alpha_ = md_gcv.best_params_['alpha']	# 最优参数
      print('训练完毕，最优参数为：', alpha_)
  
      r2_cv = np.zeros((len(alphas_), n_sp))	# 初始化矩阵，存放所有参数训练模型的r2
      for i in range(n_sp):
          name = 'split{}_test_score'.format(i)	# 遍历每次训练中的所有验证精度
          r2_cv[:, i] = md_gcv.cv_results_[name]
  
      plt.figure()
      plt.semilogx(alphas_, r2_cv, ':')	# 画出所有参数与r2的关系图
      plt.plot(alphas_, r2_cv.mean(axis=-1), 'k',
               label='Average across the folds', linewidth=2)	# 画出每个参数在所有模型的平均r2
      plt.axvline(alpha_, linestyle='--', color='k',
                  label='alpha: CV estimate')		# 画出最优参数所处的r2
      plt.legend()
      plt.xlabel(r'$\alpha$')
      plt.ylabel('Mean test score')
      plt.title('Mean test score on each fold: {}'.format(estimator.__name__))
      plt.axis('tight')
      plt.show()
      return md_gcv_res
  
  
  # 定义交叉验证的折数
  n_splits = 10
  
  # 定义训练集与验证集划分方法
  kf = ShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=2021)
  
  # 设定参数
  my_params = {'alpha': np.logspace(-5, 2, 100)}
  
  # 画出图像
  lasso_md = plot_cv(Lasso, std_x, y_train, my_params,
                     cv=kf, n_jobs=-1, scoring='r2', n_sp=n_splits)
  print('Lasso回归的R2：', lasso_md.score(std.transform(x_test), y_test))
  ridge_md = plot_cv(Ridge, std_x, y_train, my_params,
                     cv=kf, n_jobs=-1, scoring='r2', n_sp=n_splits)
  print('岭回归的R2', ridge_md.score(std.transform(x_test), y_test))
  

• 展示的结果如下：
  
  
  

Pipeline的使用

• 为了寻找最优模型，一般会选择多个模型进行比较，而过程中有许多步骤是可以重复性进行的，因此可以考虑流水线式的操作，使用Pipeline（上面模型比较中就是使用了make_pipeline，也是其中一种方式）
• 用Pipeline和GridSearchCV进行调参
  • 注意，下面加入了StandardScaler()，实际上自动将输入的train先进行标准化，输入test时会自动转化，具体可以看sklearn.base.TransformerMixin的框架思路。
  • 上述模型（Lasso回归和岭回归）都是需要先标准化处理，常用：StandardScaler。

  # 初始化Pipe
  md = Pipeline([('std', StandardScaler()),
                 ('clf', GridSearchCV(Lasso(), my_params, cv=kf,
                                      n_jobs=-1, scoring='r2', verbose=0))])
  md.fit(x_train, y_train)	# 拟合
  md.score(x_test, y_test)	# 评价
  
  # 访问模型的最优参数
  print(md['clf'].best_params_)
  

九、本文内容的所有代码

import numpy as np
import pandas as pd
from sklearn.datasets import load_boston
import os
from sklearn.linear_model import LinearRegression, Lasso, Ridge
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline, Pipeline
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.metrics import r2_score
from sklearn.model_selection import GridSearchCV, ShuffleSplit
from sklearn.linear_model import LassoCV, RidgeCV
import matplotlib.pyplot as plt


# 定义向前逐步回归
class ForwardPipeline(BaseEstimator, RegressorMixin):
    def __init__(self):
        """
        初始化，参考BaseEstimator的定义方式
        """
        # self.feature = None
        self.md = None
        self.transform = None
        self.best_feature = None
        self.mean = 0
        self.std = 0

    def fit(self, x_, y_):
        """
        拟合最优模型
        """
        # 记录变量标签名
        # self.feature = list(x_.keys())

        # 数据标准化
        # self.transform = StandardScaler()
        # self.transform.fit(x_)
        # self.mean = np.mean(x_, axis=0)
        # self.std = np.std(x_, axis=0)
        # std_x = (x_ - self.mean) / self.std
        # std_x = pd.DataFrame(std_x, columns=self.feature)
        std_x = x_

        # 求最优参数
        self.best_feature = self.forward_reg(std_x, y_)
        self.md = sm.OLS(y_, std_x[self.best_feature]).fit()
        return self

    def predict(self, x_):
        # std_x = pd.DataFrame(self.transform.transform(x_), columns=self.feature)
        # std_x = (x_ - self.mean) / self.std
        std_x = x_
        y_pred = self.md.predict(std_x[self.best_feature])
        return y_pred

    def score(self, x_, y_, sample_weight=None):
        """
        计算R2，参考RegressorMixin的定义方式，对RegressorMixin的方法进行重写
        """
        y_pred = self.predict(x_)
        return r2_score(y_, y_pred, sample_weight=sample_weight)

    @staticmethod
    def forward_reg(x_, y_, cons=True):
        md_scores = list()
        choose = list()
        train_x = x_
        if cons:
            train_x = sm.add_constant(x_, prepend=True)
            choose.append('const')
        feature = list(x_.keys())
        md_scores.append(sm.OLS(y_, np.ones(len(y_))).fit().aic)  # 模型不含特征变量时
        loop = 0
        while feature:
            best_score = float('inf')  # 初始化为无穷大
            choose_idx = 0  # 初始化单轮选择的变量名
            for i in range(len(feature)):
                score = sm.OLS(y_, train_x[choose + [feature[i]]]).fit().aic
                if score < best_score:
                    best_score = score
                    choose_idx = i
            choose.append(feature.pop(choose_idx))
            md_scores.append(best_score)
            loop += 1
            # print('第{}轮，加入变量{}，AIC为：{}'.format(loop, choose[-1], best_score))
        best_score = float('inf')
        best_idx = 0
        for idx, md_score in enumerate(md_scores):
            if md_score < best_score:
                best_idx = idx
                best_score = md_score
        # print('最终选取特征变量为：{}'.format(choose[1:best_idx + 1]))
        return choose[1:best_idx + 1]


# 定义可视化函数，参数与精度关系的可视化（基于GridSearchCV）
def plot_cv(estimator, x_, y_, params_, cv, n_jobs=-1, scoring='r2', n_sp=10):
    md = estimator()
    md_gcv = GridSearchCV(md, params_, cv=cv, n_jobs=n_jobs, scoring=scoring, verbose=0)
    md_gcv_res = md_gcv.fit(x_, y_)
    alphas_ = params_['alpha']
    alpha_ = md_gcv.best_params_['alpha']
    print('训练完毕，最优参数为：', alpha_)

    mse_cv = np.zeros((len(alphas_), n_sp))
    for i in range(n_sp):
        name = 'split{}_test_score'.format(i)
        mse_cv[:, i] = md_gcv.cv_results_[name]

    plt.figure()
    plt.semilogx(alphas_, mse_cv, ':')
    plt.plot(alphas_, mse_cv.mean(axis=-1), 'k',
             label='Average across the folds', linewidth=2)
    plt.axvline(alpha_, linestyle='--', color='k',
                label='alpha: CV estimate')
    plt.legend()
    plt.xlabel(r'$\alpha$')
    plt.ylabel('Mean test score')
    plt.title('Mean test score on each fold: {}'.format(estimator.__name__))
    plt.axis('tight')
    plt.show()
    return md_gcv_res


if __name__ == '__main__':
    # 设定项目路径
    os.chdir('D:/NewProject/datawhale/ensemble/task3')
    np.set_printoptions(suppress=True)
    np.random.seed(2021)

    # 加载数据集
    data = load_boston()
    print(data.keys())  # 查看数据集的字段

    # 选定变量
    x = data['data']  # 返回的是ndarray类型
    y = data['target']

    # 转化成dataframe结构
    d_x = pd.DataFrame(x, columns=data['feature_names'])
    d_y = pd.DataFrame(y, columns=['PRICE'])
    data_set = pd.concat([d_x, d_y], axis=1)

    # 分割数据集
    x_train, x_test, y_train, y_test = train_test_split(d_x, d_y, test_size=0.2)

    # 普通ols线性回归模型
    md1 = make_pipeline(LinearRegression())
    md1.fit(x_train, y_train)
    print('普通ols回归的R2：', md1.score(x_test, y_test))

    # 岭回归，超参数alpha（又称为岭迹）
    md2 = make_pipeline(StandardScaler(), Ridge(alpha=0.5))
    md2.fit(x_train, y_train)
    print('岭回归的R2：', md2.score(x_test, y_test))

    # Lasso回归，超参数alpha（正则化系数）
    md3 = make_pipeline(StandardScaler(), Lasso(alpha=0.5))
    md3.fit(x_train, y_train)
    print('Lasso回归的R2：', md3.score(x_test, y_test))

    # PCA降维后回归，这里降维后选取的主成分个数也算是超参数
    md4 = make_pipeline(StandardScaler(), PCA(n_components=0.8), LinearRegression())
    md4.fit(x_train, y_train)
    print('PCA降维后回归的R2：', md4.score(x_test, y_test))

    # 向前逐步回归，这里并未带有超参数
    md5 = make_pipeline(ForwardPipeline())
    md5.fit(x_train, y_train)
    print('向前逐步回归的R2：', md5.score(x_test, y_test))

    # 自定义分割训练集和验证集的方法
    kf = ShuffleSplit(n_splits=10, test_size=0.2, random_state=2021)

    # 对Lasso回归进行调参
    lasso = Lasso()
    lasso_param = {'alpha': [i for i in np.arange(0.1, 1, 0.01)]}
    lasso_gcv = GridSearchCV(lasso, lasso_param, n_jobs=-1, cv=kf, scoring='r2', verbose=1)
    lasso_gcv_res = lasso_gcv.fit(x_train, y_train)
    print('Lasso回归最优参数：', lasso_gcv_res.best_params_)
    print('Lasso回归得分：', lasso_gcv_res.score(x_test, y_test))

    # 对岭回归进行调参
    ridge = Ridge()
    ridge_param = {'alpha': [i for i in np.arange(0.1, 1, 0.05)]}
    ridge_gcv = GridSearchCV(ridge, ridge_param, n_jobs=-1, cv=kf, scoring='r2', verbose=1)
    ridge_gcv_res = ridge_gcv.fit(x_train, y_train)
    print('岭回归回归最优参数：', ridge_gcv_res.best_params_)
    print('岭回归得分：', ridge_gcv_res.score(x_test, y_test))

    # 比较模型
    m1 = LinearRegression()
    m1.fit(x_train, y_train)

    m2 = Lasso(alpha=lasso_gcv.best_params_['alpha'])
    m2.fit(x_train, y_train)

    m3 = Ridge(alpha=ridge_gcv.best_params_['alpha'])
    m3.fit(x_train, y_train)

    print('模型结果如下：\nols回归：{}\nLasso回归：{}\n岭回归：{}'
          .format(m1.score(x_test, y_test), m2.score(x_test, y_test), m3.score(x_test, y_test)))

    # 标准化处理后
    m1 = LinearRegression()
    m1_pipe = Pipeline([('sc', StandardScaler()),
                        ('md', m1)])
    m1_pipe.fit(x_train, y_train)

    m2 = Lasso()
    m2_pipe = Pipeline([('sc', StandardScaler()),
                        ('gcv', GridSearchCV(m2, lasso_param, n_jobs=-1, cv=kf, scoring='r2'))])
    m2_pipe.fit(x_train, y_train)

    m3 = Ridge()
    m3_pipe = Pipeline([('sc', StandardScaler()),
                        ('gcv', GridSearchCV(m3, ridge_param, n_jobs=-1, cv=kf, scoring='r2'))])
    m3_pipe.fit(x_train, y_train)

    print('模型结果如下：\nols回归：{}\nLasso回归：{}\n岭回归：{}'
          .format(m1_pipe.score(x_test, y_test), m2_pipe.score(x_test, y_test), m3_pipe.score(x_test, y_test)))

    # 使用LassoCV和RidgeCV调参
    std = StandardScaler()
    std_x = std.fit_transform(x_train)  # 标准化处理
    ls_ = LassoCV(eps=0.0001, n_alphas=100, cv=kf).fit(std_x, y_train)
    ridge_ = RidgeCV(alphas=np.logspace(-4, 4, 100), cv=kf).fit(std_x, y_train)
    print('Lasso回归的最优参数', ls_.alpha_)
    print('岭回归的最优参数', ridge_.alpha_)

    # 官方文档中参数与精度关系的可视化设计（基于LassoCV）
    epsilon = 1e-4  # 用于避免semilogx中x为0
    plt.figure()
    plt.semilogx(ls_.alphas_ + epsilon, ls_.mse_path_, ':')
    plt.plot(ls_.alphas_ + epsilon, ls_.mse_path_.mean(axis=-1), 'k',
             label='Average across the folds', linewidth=2)
    plt.axvline(ls_.alpha_ + epsilon, linestyle='--', color='k',
                label='alpha: CV estimate')
    plt.legend()
    plt.xlabel(r'$\alpha$')
    plt.ylabel('Mean square error')
    plt.title('Mean square error on each fold: coordinate descent ')
    plt.axis('tight')
    plt.show()

    # 定义交叉验证的折数
    n_splits = 10

    # 定义训练集与验证集划分方法
    kf = ShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=2021)

    # 设定参数
    my_params = {'alpha': np.logspace(-5, 2, 100)}

    # 画出图像
    lasso_md = plot_cv(Lasso, std_x, y_train, my_params,
                       cv=kf, n_jobs=-1, scoring='r2', n_sp=n_splits)
    print('Lasso回归的R2：', lasso_md.score(std.transform(x_test), y_test))
    ridge_md = plot_cv(Ridge, std_x, y_train, my_params,
                       cv=kf, n_jobs=-1, scoring='r2', n_sp=n_splits)
    print('岭回归的R2', ridge_md.score(std.transform(x_test), y_test))

    # Pipeline的使用
    mm3 = Pipeline([('std', StandardScaler()),
                    ('clf', GridSearchCV(Lasso(), my_params, cv=kf,
                                         n_jobs=-1, scoring='r2', verbose=0))])
    mm3.fit(x_train, y_train)
    mm3.score(x_test, y_test)

十、参考资料

1. https://github.com/datawhalechina/ensemble-learning
2. https://www.bilibili.com/video/BV1Mb4y1o7ck?t=470
3. https://blog.csdn.net/ghj786110/article/details/106711047/
4. https://zhuanlan.zhihu.com/p/71961236