ensemble learning 学习分享：回归问题

本文是学习开源项目‘ensemble learning’之后的总结。平台使用jupyter notebook。

---

项目链接：https://github.com/datawhalechina/ensemble-learning

章节学习视频链接：https://www.bilibili.com/video/BV1Mb4y1o7ck?p=1 前10p

---

机器学习的任务可分为有监督学习和无监督学习，典型的有监督学习如回归和分类，本节仅介绍回归和分类相关的问题。

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

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

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

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

用sklearn构建完整的回归项目

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

选用scikit-learn自带的数据集：Boston房价数据集，以避免在处理数据上花费较大精力，达到快速上手sklearn掌握相关操作的效果。

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

各个特征的相关解释：

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

选择度量模型性能的指标

• 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|$，平均绝对误差是一个与绝对误差损失（或L1范数损失）期望值的相应的风险指标。

• $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^​}​，解释方差得分最好的是1，值越低越差。

  以上， $\widehat{y_i}$ 是第 $i$ 个样本的预测值，$y_i$ 是相应的真实值； $\hat{y}$ 是估计的 target output，$y$ 是相应的(correct) target output。

各指标含义及代码实现具体可参见：Regression metrics

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

线性回归模型：

回归的具体概念可了解一下高尔顿关于父子身高的研究，并理解一下”向均值回归“的含义即可。

假设影响因变量（dependent variable） $Y$ 的自变量（independent variable）有 $p$ 个：$X_1,X_2,\cdot \cdot \cdot ,X_p$ , 并且设定它们之间有如下线性关系：
$$Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdot \cdot \cdot +{\beta }_p{X}_p+e,$$
其中 $e$ 是误差项， ${\beta }_0,{\beta }_1,\cdots ,{\beta }_p$ 是待估计的未知参数。自变量也称为预测变量（predictor variables）或控制变量（control variables）。

假设我们获得了响应变量 $Y$ 和预测变量 $X_1,X_2,\cdot \cdot \cdot ,X_p$ 的 $n$ 组观测值$(x_{i1},\cdots ,x_{ip},y_i)，i=1,2,\cdots ,n$，它们满足
$$y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_p{x}_{ip}+e_i.$$
用矩阵形式表示为：
$$\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_i\end{array}\right]=\left[\begin{array}{ccccc}1& x_{11}& x_{12}& \cdots & x_{1p}\\ 1& x_{21}& x_{22}& \cdots & x_{2p}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& x_{n1}& x_{n2}& \cdots & x_{np}\end{array}\right]\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\\ ⋮\\ {\beta }_i\end{array}\right]+\left[\begin{array}{c}{e}_1\\ e_2\\ ⋮\\ e_i\end{array}\right]$$
即$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{e}$。

模型的基本假定：

(1). 解释变量$X_1,X_2,\cdot \cdot \cdot ,X_p$ 是确定性变量，不是随机变量。且要求设计矩阵$\mathbf{X}$ 的秩为p+1，即满秩矩阵

(2). 随机误差项具有0均值和等方差，即
$$\{\begin{array}{lc}E(e_i)=0,& i=1,2,\cdots ,n\\ cov(e_i,e_j)=\{\begin{array}{lc}{e}^2,& i=j\\ 0,& i\ne j\end{array}& i,j=1,2,\cdots ,n\end{array}$$
​ 这个假定又被称为高斯-马尔科夫（Gauss-Markov）条件。

(3). 正态分布的假定条件：
$$\{\begin{array}{l}{e}_i\sim N(0,{\sigma }^2),i=1,2,\cdots ,n\\ e_1,e_2,\cdots ,e_n相互独立\end{array}$$
即$\mathbf{e}\sim N(0,{\sigma }^2{\text{I}}_{\text{n}})$ 。

则线性回归模型中最小二乘法表述如下：

寻找参数${\beta }_0,{\beta }_1,{\beta }_2,\cdots ,{\beta }_p$ 的估计值 ${\hat{\beta }}_0,{\hat{\beta }}_1,{\hat{\beta }}_2,\cdots ,{\hat{\beta }}_p$ 使离差平方和$Q({\beta }_0,{\beta }_1,{\beta }_2,\cdots ,{\beta }_p)=\sum _{i=1}^n(y_i-{\beta }_0-{\beta }_1{x}_{i1}-{\beta }_2{x}_{i2}-\cdots -{\beta }_p{x}_{ip}{)}^2$ 达到极小，写成矩阵形式为：$Q(\beta )=\sum _{i=1}^N||{\beta }^T{x}_i-y_i|{|}_2^2=\sum _{i=1}^N({\beta }^T{x}_i-y_i{)}^2=(Y-\beta X{)}^T(Y-\beta X)={\mathbf{e}}^{\mathbf{T}}\mathbf{e}$ ，即为损失函数。

对${\beta }_i$ 求偏导并令其为0，整理得：${\mathbf{X}}^T(\mathbf{y}-\mathbf{X}\hat{\mathbf{\beta }})=0$ ，当$(\mathbf{X}{\mathbf{X}}^{\mathbf{T}}{)}^{-1}$ 存在时，得参数的最小二乘估计 $\beta =((\mathbf{X}{\mathbf{X}}^{\mathbf{T}}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{y})$

涉及到的矩阵求导知识具体参见：矩阵求导术（上） 、矩阵求导术（下）

该部分内容参考：何晓群、刘文卿《应用回归分析》（第5版） 第3章、林建忠《回归分析与线性统计模型》 第2章

关于“在线性回归模型中，极大似然估计与最小二乘估计有什么联系与区别？”，可参见：最大似然估计和最小二乘估计的区别与联系

另外，关于最小二乘的几何解释在林建忠《回归分析与线性统计模型》 第2章中也可找到，不明白者可翻阅查看。

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
'''

树模型与线性模型的区别 决策树分类和逻辑回归分类的区别 【总结】 回归决策树和线性回归代码对比

关于KKT条件：

具体参见：如何理解拉格朗日乘子法和KKT条件？

用线性回归模型拟合数据集，要求只使用numpy

注：这里仅就一元线性回归问题做代码实现

# 一元线性回归问题
# y = w * x + b，w与b是参数

import numpy as np

# 定义损失函数
def loss_function(w, b, points):
    loss = 0
    for i in range(0, len(points)):
        x = points[i, 0]                  # x值
        y = points[i, 1]                  # 相应的y值
        loss += (y - (w * x + b)) ** 2   # MSE均方误差*样本量
        loss = loss / float(len(points))  # 损失函数—MSE均方误差
        return loss

# 实现一次对w 和 b 的更新
def step_gradient(w_current, b_current, points, learning_rate):
    w_gradient = 0
    b_gradient = 0
    N = float(len(points))            # 样本量
    for i in range(0, len(points)):
        x = points[i, 0]
        y = points[i, 1]
        w_gradient += -(2 / N) * (y - (w_current * x + b_current)) * x  # 计算w的梯度，即对w求导，最小二乘
        b_gradient += -(2 / N) * (y - (w_current * x + b_current))      # 计算b的梯度，即对b求导，最小二乘
        new_w = w_current - (learning_rate * w_gradient)                 # 最小二乘+一次梯度下降法得到的w的估计值
        new_b = b_current - (learning_rate * b_gradient)                 # 最小二乘+一次梯度下降法得到的b的估计值
    return [new_w, new_b]


# 调用step_gradient()，对w,b循环更新num_iterations次，返回w,b
def gradient_descent_runner(points, starting_w, starting_b, learning_rate, num_iterations):
    w = starting_w
    b = starting_b
    for i in range(num_iterations):
        w, b = step_gradient(w, b, np.array(points), learning_rate)
    return [w, b]


def main():
    points = np.genfromtxt("data.csv", delimiter=",")   # 读入csv文件,成ndarray
    learning_rate = 0.0001     # 初始化学习率
    initial_w = 0              # 初始化带估参数w
    initial_b = 0              # 初始化待估参数b
    num_iterations = 1000      # 设置迭代次数
    print("Starting gradient descent at w = {0},b = {1},loss = {2}".format(initial_w, initial_b,
                                                                           loss_function(initial_w, initial_b, points)))
    print("Running...")
    w, b = gradient_descent_runner(points, initial_w, initial_b, learning_rate, num_iterations)
    print("After {0} iterations w = {1},b = {2},loss = {3}".format(num_iterations, w, b, loss_function(w, b, points)))


if __name__ == '__main__':
    main()
    
'''
Starting gradient descent at w = 0,b = 0,loss = 10.053342197544191
Running...
After 1000 iterations w = 1.4777440851894448,b = 0.08893651993741344,loss = 2.693563420287152
'''

多元回归的代码实现具体参见：纯python（numpy）实现多元线性回归 及 线性回归理解（附纯python实现）