XGBoost -Regression

mzx8181

已于 2023-08-07 19:31:48 修改

阅读量232

点赞数

文章标签： python

于 2023-08-01 16:42:12 首次发布

本文链接：https://blog.csdn.net/mzx8181/article/details/132043407

版权

1 Classification and Regression Tree (CART)：

CART是一种决策树模型，可以用于分类或者回归.CART学习算法共分为两步：第一步是决策树生成：用训练数据生成决策树，生成的树尽可能大；第二步决策树剪枝：基于损失函数最小化进行剪枝，用验证数据对生成的数据剪枝。针对回归树，使用平方损失来表示回归树对于训练数据的预测误差，平方误差最小处即为最优输出值

$\sum\limits_{x_i \in R_m}(y_i - f(x_i))^2$

回归树采用采用启发式的方法，根据特征分裂成左子树和右子树。如，对于第 $j$ 个变量 $x_j$ 作为划分特征， $x_j=s$ 为其分裂点，那么样本就被划分为：

$R_1(j,s) = \{x|x_i\leqslant {s}\}$
$R_2(j,s) = \{x|x_i>s\}$
具体求解最优划分特征 $j$ 和最优划分点 $s$

$\min\limits_{j,s}[\min\limits_{c_1}\sum\limits_{x_i \in R_1}(y_i - c_1)^2+ \sum\limits_{x_i \in R_2}(y_i - c_2)^2]$

只需要求解所有特征的所有切分点，就可以找到最优切分特征和切分点，得到最终的回归树

2 Extreme Gradient Boosting (XGBoost)

（1）XGBoost算法思想
XGBoost是Chen 和 Guestrin 提出的一种集成学习模型，其借助了 CART 回归树的思想，并在其基础上进行了改进。算法采用Gradient Boosting的思想，每一个基学习器重点关注前一个基学习器不足的地方进行训练，串行的训练多个模型逐步逼近降低损失，最终根据样本特征将每棵树落到对应的叶子结点上，将每个叶子节点上的分数相加，即为预测值。
Boosting是串行地训练多个模型，逐步逼近降低损失

$\hat{y}=\sum^K_{k=1}f_k(x_i)$

（2）XGBoost算法原理
XGBoost目标函数定义为
$Obj=\sum^n_{i=1}l(y_i,\hat{y_i})+\sum^K_{k=1}\Omega(f_k)$

由于XGBoost是在前k-1棵树的基础上进行预测的，因此前k-1棵树的预测结果和模型复杂度是已知的，只有第k棵树的预测值是未知的，因此，生成第k棵树之后，预测值就可以表示为

$\hat{y_i}^{(t)}=\hat{y_i}^{(t-1)}+f_t(x_i)$

由于 $\sum_{1}^{k-1}\Omega(f_k)$ 可以视为常数省略，因此，目标函数就可以表示为

$Obj=\sum^n_{i=1}l(y_i,\hat{y_i}^{(t-1)}+f_t(x_i))+\Omega(f_k)$
利用泰勒二阶级数展开目标函数，可以得到目标函数为

$\approx \sum_{i=1}^{n}[l(y_i,\hat{y_i}^{t-1})+\partial_{\hat{y_i}^{t-1}}l(y_i,\hat{y_i}^{t-1})f_t(x_i)+\frac{1}{2}\partial^{2}_{\hat{y_i}^{t-1}}l(y_i,\hat{y_i}^{t-1})f^2_t(x_i)]+\Omega(f_t)$

令：
$g_i = \partial_{\hat{y_i}^{t-1}}l(y_i,\hat{y_i}^{t-1})$

$h_i=\partial^{2}_{\hat{y_i}^{t-1}}l(y_i,\hat{y_i}^{t-1})$
目标函数化简为：

$\approx [g_if_k(x_i)+\frac{1}{2}h_if^2_k(x_i)]+\Omega(f_k)$
为了将目标函数与树的结构相结合，
将 $\omega_{q(x)}$ 定义为叶子结点 $q$ 的分数
样本 $x_i$ 落在相应叶子结点上的预测值为 $f_t(x_i) = w_{q(x_i)}$
将第k棵树的模型复杂度与树结构相结合：

$\Omega(f_k)=\gamma{T}+\frac{1}{2}\lambda\sum_{j=1}^{T}w^2_j$

将目标函数Obj用树结构进行表示：

$Obj=\sum_{i=1}^n[g_iw_{q(x_i)}+\frac{1}{2}h_iw_{q(x_i)}^2]+\gamma{T}+\frac{1}{2}\lambda\sum_{j=1}^Tw_j^2$
其中，
$T$ 表示为叶子结点的个数， $\gamma$ 和 $\lambda$ 是超参数，用于控制模型复杂度

令

$G_j=\sum_{i\epsilon I_j}g_j$

$H_j=\sum_{i\epsilon I_j}h_j$

目标函数化简为：

$\approx \sum_{j=1}^T[G_jw_i+\frac{1}{2}(H_j+\lambda)w_i^2]+\gamma{T}$

假设已知树的结构q，那么可以通过求解关于 $\omega$ 的一元二次函数的最小值，求出最优的 $\omega$ 及对应的最大增益，即

$w^*_j=-\frac{G_J}{2 \times\frac{1}{2} (H_j+\lambda)}=-\frac{G_j}{H_j+\lambda}$

代入目标函数，

$Obj=-\frac{1}{2}\sum_{j=1}^T\frac{G^2_j}{H_j+\lambda}+\gamma{T}$ (16)

采用贪心算法的思路，对当前树中的每一个叶子节点，用不同特征尝试进行分割，并用下面的函数计算分裂前和分裂后的增益分数，找到最优的分割方案

$Gain=Obj_{L+R}-(Obj_L+Obj_R)$

$=[-\frac{1}{2}\frac{(G_L+G_R)^2}{H_L+H_R+\lambda}+\gamma{T}]-[-\frac{1}{2}(\frac{G_L^2}{H_L+\lambda}+\frac{G_R^2}{H_R+\lambda})+\gamma{(T+1)}]$

$=\frac{1}{2}[\frac{G_L^2}{H_L+\lambda}+\frac{G_R^2}{H_R+\lambda}-\frac{G_L+G_R)^2}{H_L+H_R+\lambda}]-\gamma$

3 XGBoost Source Code

import pandas as pd
import numpy as np


class XGB:

    def __init__(self,
                 base_score=0.5,
                 max_depth=3,
                 n_estimators=10,
                 learning_rate=0.1,
                 reg_lambda=1,
                 gamma=0,
                 min_child_sample=None,
                 min_child_weight=1,
                 objective='linear'):

        self.base_score = base_score  # 最开始时给叶子节点权重所赋的值，默认0.5，迭代次数够多的话，结果对这个初值不敏感
        self.max_depth = max_depth  # 最大树深度
        self.n_estimators = n_estimators  # 树的个数
        self.learning_rate = learning_rate  # 学习率，别和梯度下降里的学习率搞混了，这里是每棵树要乘以的权重系数
        self.reg_lambda = reg_lambda  # L2正则项的权重系数
        self.gamma = gamma  # 正则项中，叶子节点数T的权重系数
        self.min_child_sample = min_child_sample  # 每个叶子节点的样本数
        self.min_child_weight = min_child_weight  # 每个叶子节点的Hessian矩阵和，下面代码会细讲
        self.objective = objective  # 目标函数，可选linear和logistic
        self.tree_structure = {}  # 用一个字典来存储每一颗树的树结构

    def xgb_cart_tree(self, X, w, m_dpth):
        '''
        递归的方式构造XGB中的Cart树
        X：训练数据集
        w：每个样本的权重值，递归赋值
        m_dpth：树的深度
        '''

        # 边界条件：递归到指定最大深度后，跳出
        if m_dpth > self.max_depth:
            return

        best_var, best_cut = None, None

        # 这里增益的初值一定要设置为0，相当于对树做剪枝，即如果算出的增益小于0则不做分裂
        max_gain = 0
        G_left_best, G_right_best, H_left_best, H_right_best = 0, 0, 0, 0

        # 遍历每个变量的每个切点，寻找分裂增益gain最大的切点并记录下来
        # 第一层循环是提取每个特征的列名item
        for item in [x for x in X.columns if x not in ['g', 'h', 'y']]:
            # list(set(X[item])) 提取每个特征X[item]下数据的唯一值，并构建一个去重数据的列表
            # cut是遍历列表中的每个值
            for cut in list(set(X[item])):

                # 考虑叶子节点最小样本数：如果指定了min_child_sample则限制分裂后叶子节点的样本数都不能小于指定值
                if self.min_child_sample:
                    # 对X数据集的X[item]列进行条件筛选
                    # 如果在cut切分点下有任意一遍的叶子节点样本数小于min_child_sample，就跳出本层循环，即当前遍历到的cut值，进行下一个cut值的条件判断
                    # \是为了增强代码可读性而进行换行的符号
                    if (X.loc[X[item] < cut].shape[0] < self.min_child_sample) \
                            | (X.loc[X[item] >= cut].shape[0] < self.min_child_sample):
                        continue

                # 选取X中满足条件X[item] < cut的行，并且选择了这些行中的 'g' 列（loss一阶导数）
                # .sum() 方法对这些选定的 'g' 列进行求和操作，得到结果 G_left（loss一阶导数的和）
                G_left = X.loc[X[item] < cut, 'g'].sum()
                G_right = X.loc[X[item] >= cut, 'g'].sum()
                # 选取X中满足条件X[item] < cut的行，并且选择了这些行中的 'h' 列（loss二阶导数）
                # .sum() 方法对这些选定的 'g' 列进行求和操作，得到结果 G_left（loss二阶导数的和）
                H_left = X.loc[X[item] < cut, 'h'].sum()
                H_right = X.loc[X[item] >= cut, 'h'].sum()

                # min_child_weight 每个叶子节点的Hessian矩阵和，用来控制纯度，防止过拟合
                # min_child_weight在这里起作用，指的是每个叶子节点上的H，即目标函数二阶导的加和
                # 当目标函数为linear，即1/2*(y-y_hat)**2时，它的二阶导是1，那min_child_weight就等价于min_child_sample
                # 当目标函数为logistic，其二阶导为sigmoid(y_hat)*(1-sigmoid(y_hat))，可理解为叶子节点的纯度
                if self.min_child_weight:
                    if (H_left < self.min_child_weight) | (H_right < self.min_child_weight):
                        continue

                # 增益gain的计算，找到每个特征下的gain值最大值的分裂点
                gain = G_left ** 2 / (H_left + self.reg_lambda) + \
                       G_right ** 2 / (H_right + self.reg_lambda) - \
                       (G_left + G_right) ** 2 / (H_left + H_right + self.reg_lambda)
                gain = gain / 2 - self.gamma
                # 函数定义了初始的max_gain=0
                if gain > max_gain:
                    # best_var 最佳分割变量
                    # best_cut 最佳分裂点
                    best_var, best_cut = item, cut
                    max_gain = gain
                    G_left_best, G_right_best, H_left_best, H_right_best = G_left, G_right, H_left, H_right

        # 如果遍历完找不到可分列的点，则返回None
        if best_var is None:
            return None

        # 给每个叶子节点上的样本分别赋上相应的权重值
        # 将数据集X按照最佳分割变量best_var的最佳分裂点best_cut分成左子节点和右子节点

        # 左子节点
        # 筛选出左子节点行索引，并通过index.tolist()转化为列表赋给id_left
        id_left = X.loc[X[best_var] < best_cut].index.tolist()
        # 计算左子节点的权重
        w_left = - G_left_best / (H_left_best + self.reg_lambda)

        # 右子节点
        # 筛选出右子节点行索引，并通过index.tolist()转化为列表赋给id_right
        id_right = X.loc[X[best_var] >= best_cut].index.tolist()
        # 计算左子节点的权重
        # 目的是使得梯度（G_right_best）和Hessian矩阵和（H_right_best）越大，权重越小，从而降低损失函数的值
        w_right = - G_right_best / (H_right_best + self.reg_lambda)

        # 将左子节点和右子节点的权重w_left和w_right分别赋值给权重数组w中相应的索引位置，更新权重数组
        w[id_left] = w_left
        w[id_right] = w_right

        # 用俄罗斯套娃式的json串把树的结构给存下来
        # best_var, best_cut当前节点的划分特征和划分点，是字典的顶层键
        tree_structure = {(best_var, best_cut): {}}

        # ('left', w_left)是左子树的键 w_left是左子树的权重
        # 递归调用 xgb_cart_tree 方法构建左子树的结构
        # X.loc[id_left]是基于左子树样本的子数据集
        # w是权重数组
        # m_dpth+1当前树的深度加1，指的是左子树的深度
        tree_structure[(best_var, best_cut)][('left', w_left)] = self.xgb_cart_tree(X.loc[id_left], w, m_dpth + 1)

        # ('right', w_right)是左子树的键 w_right是右子树的权重
        # 递归调用 xgb_cart_tree 方法构建右子树的结构
        # X.loc[id_right]是基于右子树样本的子数据集
        # w是权重数组
        # m_dpth+1当前树的深度加1，指的是右子树的深度
        # 左子树和右子树的深度都是在当前的深度上+1，因此相等
        tree_structure[(best_var, best_cut)][('right', w_right)] = self.xgb_cart_tree(X.loc[id_right], w, m_dpth + 1)

        return tree_structure

    # 定义一个函数计算目标函数一阶导
    # 预测值y_hat和真实值Y的差值y_hat-Y是平方损失函数的一阶导数
    def _grad(self, y_hat, Y):
        '''
        计算目标函数的一阶导
        支持linear和logistic
        y_hat 预测值
        Y 真实值
        '''

        if self.objective == 'logistic':
            # 首先将预测值 y_hat 应用逻辑函数，得到概率值 p
            # 然后计算概率值 p 与真实值 Y 的差异，即 p - Y，作为逻辑回归的一阶导数
            y_hat = 1.0 / (1.0 + np.exp(-y_hat))
            return y_hat - Y
        elif self.objective == 'linear':
            # 如果是线性回归就直接求差值作为一阶导数
            return y_hat - Y
        else:
            raise KeyError('objective must be linear or logistic!')

    # 定义一个函数计算目标函数二阶导
    def _hess(self, y_hat, Y):
        '''
        计算目标函数的二阶导
        支持linear和logistic
        '''

        if self.objective == 'logistic':
            # 转化概率值
            y_hat = 1.0 / (1.0 + np.exp(-y_hat))
            # 逻辑回归目标函数的二阶导数
            return y_hat * (1.0 - y_hat)
        elif self.objective == 'linear':
            # 线性回归的二阶导数是常数，即每个样本的二阶导数都是1
            return np.array([1] * Y.shape[0])
        else:
            raise KeyError('objective must be linear or logistic!')

    def fit(self, X: pd.DataFrame, Y):
        '''
        根据训练数据集X和Y训练出树结构和权重
        X 特征值
        Y 目标值
        '''

        # 检查X和Y是否有相同的行数
        if X.shape[0] != Y.shape[0]:
            raise ValueError('X and Y must have the same length!')

        # 重置X的索引，丢弃旧的索引，生成一个新的从0开始的索引
        X = X.reset_index(drop='True')
        # values转化为numpy数组
        Y = Y.values
        # 这里根据base_score参数设定权重初始值
        y_hat = np.array([self.base_score] * Y.shape[0])
        # 迭代，依次形成训练树
        for t in range(self.n_estimators):
            print('fitting tree {}...'.format(t + 1))

            # 一阶导数g、二阶导数h，导数值由_grad,_hess计算得到
            # y_hat预测值，Y目标变量
            X['g'] = self._grad(y_hat, Y)
            X['h'] = self._hess(y_hat, Y)

            # f_t用于存储当前树的预测结果
            f_t = pd.Series([0] * Y.shape[0])
            # 调用xgb_cart_tree根据X、f_t和索引构建树结构，并存储在self.tree_structure（字典）中
            self.tree_structure[t + 1] = self.xgb_cart_tree(X, f_t, 1)

            # 更新预测值 y_hat，通过加上学习率乘以当前树的预测结果 f_t
            y_hat = y_hat + self.learning_rate * f_t

            print('tree {} fit done!'.format(t + 1))

        print(self.tree_structure)

    def _get_tree_node_w(self, X, tree, w):
        '''
        以递归的方法，把树结构解构出来，把每个叶子结点的权重值赋到w上面
        tree是一个字典对象，表示整个树的结构，即tree_structure
        tree的每个键表示一个节点，对应的值是该节点的子树结构
        '''

        if not tree is None:
            # 提取树结构中的第一个节点，即根节点
            k = list(tree.keys())[0]
            # 获取树的划分特征和切分点
            var, cut = k[0], k[1]
            # 根据切分点将输入数据X分为左子树
            X_left = X.loc[X[var] < cut]
            # 获取左子树的索引
            id_left = X_left.index.tolist()
            # # 根据切分点将输入数据X分为右子树
            X_right = X.loc[X[var] >= cut]
            # 获取右子树的索引
            id_right = X_right.index.tolist()
            # 遍历树结构根节点的子节点的键
            for kk in tree[k].keys():
                if kk[0] == 'left':
                    # 存储左子树的树结构
                    tree_left = tree[k][kk]
                    # 将权重存在左子树对应的样本权重数组中
                    w[id_left] = kk[1]
                elif kk[0] == 'right':
                    # 存储右子树的树结构
                    tree_right = tree[k][kk]
                    # 将权重存在右子树对应的样本权重数组中
                    w[id_right] = kk[1]

            # 递归的调用该函数来遍历树的节点并赋值权重
            self._get_tree_node_w(X_left, tree_left, w)
            self._get_tree_node_w(X_right, tree_right, w)

    def predict_raw(self, X: pd.DataFrame):
        '''
        根据训练结果预测
        返回原始预测值
        '''
        # 重置索引，从0开始
        X = X.reset_index(drop='True')
        # 创建一个长度为 X 数据框的行数的 Pandas Series 对象，并初始化所有元素的值为 self.base_score
        Y = pd.Series([self.base_score] * X.shape[0])

        # 遍历每棵树
        for t in range(self.n_estimators):
            tree = self.tree_structure[t + 1]
            # 初始化一个与输入数据X行数相同长度的空的Pandas Series对象y_t，用于存储当前树的预测值
            y_t = pd.Series([0] * X.shape[0])
            # 将每个样本的预测值赋值给 y_t
            self._get_tree_node_w(X, tree, y_t)
            # 将学习率self.learning_rate乘以y_t，并将结果与Y相加，更新Y的值
            Y = Y + self.learning_rate * y_t

        return Y

    def predict_prob(self, X: pd.DataFrame):
        '''
        当指定objective为logistic时，输出概率要做一个logistic转换
        '''

        Y = self.predict_raw(X)
        sigmoid = lambda x: 1 / (1 + np.exp(-x))
        Y = Y.apply(sigmoid)

        return Y