基于XGBoost算法的集成学习

时间序列相关参考文章：
时间序列预测算法—ARIMA
时间序列预测算法—Prophet
时间序列预测算法—LSTM
长时间序列预测算法—Informer
时间序列分类任务—tsfresh
XGBoost时间序列预测
有季节效应的非平稳序列分析
python时间序列处理
时间序列异常值检测方法
时间序列异常值处理方法

一、XGBoost原理

1.1 提升方法（Boosting）

提升⽅法使⽤加法模型和前向分步算法。

• 加法模型
  

  其中，b(x;${\gamma }_m$)为基函数，${\gamma }_m$为基函数的参数，${\beta }_m$为基函数的系数。在给定训练数据及损失函数 L 的条件下，学习加法模型 𝑓(𝑥) 成为经验风险极小化问题：

  前向分步算法求解这⼀优化问题的思路：因为学习的是加法模型，可以从前向后，每⼀步只学习⼀个基函数及其系数，逐步逼近优化目标函数式，则可以简化优化复杂度。具体地，每步只需优化如下损失函数：

• 前向分步算法
  

  前向分步算法将同时求解从𝑚 = 1到 𝑀 所有参数的优化问题简化为逐次求解各个的优化问题。

1.2 提升决策树 （BDT）

  以决策树为基函数的提升方法为提升决策树。提升决策树模型可以表示为决策树的加法模型：

  提升决策树采⽤前向分步算法。⾸先确定初始提升决策树f0(x)=0，第m步的模型是:

  通过经验风险极小化确定下⼀棵决策树的参数：

  提升决策树使⽤以下前向分步算法：

  回归问题的提升决策树算法：

1.3 梯度提升决策树 （GBDT）

  梯度提升算法使⽤损失函数的负梯度在当前模型的值，作为回归问题提升决策树算法中残差的近似值，拟合⼀个回归树。

1.4 极限梯度提升（XGBoost）

  XGBoost 对GBDT进行了一系列优化，比如损失函数进行了二阶泰勒展开，目标函数加入正则化，支持并行和默认缺失值处理等，在可扩展性和训练速度上有了巨大的提升，但核心思想没有大的变化。作为Boosting算法，XGBoost中自然包含Boosting三要素:

• 损失函数L(x,y)：用以衡量模型预测结果与真实结果的差异。
• 弱评估器f(x)：(一般为)决策树，不同的Boosting算法使用不同的建树过程。
• 综合集成结果H(x)：即集成算法具体如何输出集成结果。

XGBoost也遵循Boosting算法的流程：

• 依据上一个弱评估器$f(x{)}_{k-1}$的结果，计算损失函数L，并使用L自适应地影响下一个弱评估器$f(x{)}_k$的构建。集成模型输出的结果，受到整体所有弱评估器$f(x{)}_0$~$f(x{)}_k$的影响。

1.4.1 XGBoost改进

  XGBoost在此基础上做了众多关键的改进，综合来看，这些改进都是基于XGBoost中两种非常关键的思想实现的：

• 第一，实现精确性和复杂度之间的平衡

  树的集成模型是机器学习中最为强大的学习器之一，这一族学习器的特点是精确性好、适用于各种场景，但运行缓慢、且过拟合风险很高，因此从学习单一决策树时起，就提供丰富的剪枝策略，目的就是为了降低各种树模型的模型复杂度，从而控制住过拟合。树模型的学习能力与过拟合风险之间的平衡，就是预测精确性与模型复杂度之间的平衡，也是 经验风险(损失函数)与结构风险(模型复杂度) 之间的平衡，这一平衡对决策树以及树的集成模型来说是永恒的议题。
  在过去，我们总是先建立效果优异的模型，再依赖于手动剪枝来调节树模型的复杂度，但在XGBoost中，精确性与复杂度会在训练的每一步被考虑到。主要体现在:
（1）XGBoost为损失函数L(y,$\hat{y}$)加入结构风险项，构成目标函数O(y,$\hat{y}$)
  在AdaBoost与GBDT当中，我们的目标是找到损失函数L(y,$\hat{y}$)的最小值，也就是让预测结果与真实结果差异最小，这一流程只关心精确性、不关心复杂度和过拟合情况。为应对这个问题，XGBoost从决策树的预剪枝流程、逻辑回归、岭回归、Lasso等经典算法的抗过拟合流程吸取经验，在损失函数中加入了控制过拟合的结构风险项，并将L(y,$\hat{y}$)+结构风险定义为目标函数O(y,$\hat{y}$)。
  这一变化让XGBoost在许多方面都与其他Boosting算法不同:例如，XGBoost是向着令目标函数最小化的目标进行训练，而不是令损失函数最小化的方向。再比如，XGBoost会优先利用结构风险中的参数来控制过拟合，而不像其他树的集成模型一样依赖于树结构参数(例如max depth,min impurity decrease等)
（2）使用全新不纯度衡量指标，将复杂度纳入分枝规则
  在之前学过的算法当中，无论Boosting流程如何进化，建立单棵决策树的规则基本都遵循我们曾经学过的CART树流程，在分类树中，我们使用信息增益(information gain)来衡量叶子的质量，在回归树中，我们使用MSE或者RMSE来衡量叶子的质量。这一流程有成熟的剪枝机制、预测精度高、能够适应各种场景，但却可能建立复杂度很高的树。
  为实现精确性与复杂度之间的平衡，XGBoost重新设定了分枝指标【结构分数】(原论文中Structure Score，也被称为质量分数Quality Score)，以及基于结构分数的【结构分数增益】(Gain of structure score)，结构分数增益可以逼迫决策树向整体结构更简单的方向生长。这一变化让XGBoost使用与传统CART略有区别的建树流程，同时在建树过程中大量使用残差(Residuals)或类残差对象作为中间变量，因此XGBoost的数学过程比其他Boosting算法更复杂。

• 第二，极大程度地降低模型复杂度、提升模型运行效率

  在任意决策树的建树过程中，都需要对每一个特征上所有潜在的分枝节点进行不纯度计算，当数据量巨大时这一计算将消耗巨量的时间，因此树集成模型的关键缺点之一就是计算缓慢，而这一缺点在实际工业环境当中是相当致命的。为了提升树模型的运算速度、同时又不极大地伤害模型的精确性，XGBoost使用多种优化技巧来实现效率提升:
1.使用估计贪婪算法、平行学习、分位数草图算法等方法构建了适用于大数据的全新建树流程
2.使用感知缓存访问技术与核外计算技术，提升算法在硬件上的运算性能
3.引入Dropout技术，为整体建树流程增加更多随机性、让算法适应更大数据
除此之外，XGBoost还保留了部分与梯度提升树类似的属性，包括：

• 弱评估器的输出类型与集成算法输出类型不一致

  对于AdaBoost或随机森林算法来说，当集成算法执行的是回归任务时，弱评估器也是回归器当集成算法执行分类任务时，弱评估器也是分类器。但对于GBDT以及基于GBDT的复杂Boosting算法们而言，无论集成算法整体在执行回归/分类/排序任务，弱评估器一定是回归器。GBDT通过sigmoid或softmax函数输出具体的分类结果，但实际弱评估器一定是回归器XGBoost也是如此。

• 拟合负梯度，且当损失函数是0.5倍MSE时，拟合残差

  任意Boosting算法都有自适应调整弱评估器的步骤。在GBDT当中，每次用于建立弱评估器的是样本X以及当下集成输出H(xi)与真实标签y之间的伪残差(也就是负梯度)。当损失函数是MSE时，负梯度在数学上等同于残差(Residual)，因此GBDT是通过拟合残差来影响后续弱评估器结构。XGBoost也是依赖于拟合残差来影响后续弱评估器结构。

• 抽样思想

  GBDT借鉴了大量Bagging算法中的抽样思想，XGBoost也继承了这一属性，因此在XGBoost当中，我们也可以对样本和特征进行抽样来增大弱评估器之间的独立性。

1.4.2 XGBoostcsklearn实现

#!pip install xgboost #安装xgboost库
import xgboost as xgb 
from xgboost import XGBRegressor
from sklearn.model_selection import cross_validate, KFold
from sklearn.model_selection import train_test_split
import pandas as pd

data = pd.read_csv("train_encode.csv",index_col=0)
data.head()

X = data.iloc[:,:-1]
y = data.iloc[:,-1]
#sklearn普通训练代码三步走：实例化，fit，score
x_train,x_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=12)
xgb_sk = XGBRegressor(random_state=12) #实例化模型
xgb_sk.fit(x_train,y_train)
xgb_sk.score(x_test,y_test) #默认指标R2  #输出0.8458883762359619

#定义所需的交叉验证方式
cv = KFold(n_splits=10,shuffle=True,random_state=12)
result_xgb_sk = cross_validate(xgb_sk,X,y,cv=cv
                               ,scoring="neg_root_mean_squared_error" #负根均方误差
                               ,return_train_score=True
                               ,verbose=True
                               ,n_jobs=-1)
result_xgb_sk

def RMSE(result,name):
    return abs(result[name].mean())
train_score=RMSE(result_xgb_sk,"train_score")
test_score=RMSE(result_xgb_sk,"test_score")
print(f"train_score: {train_score:.2f}\ntest_score: {test_score:.2f}")
#train_score: 1158.96
#test_score: 27243.69

可以看到，在默认参数下，xgboost模型不稳定，在训练集上的RMSE为1158.96，测试集为27243.69。这说明XGBoost的学习能力的确强劲，但过拟合的情况非常严重。对XGBoost的参数略微进行调整，例如将最可能影响模型的参数之一：max_depth设置为一个较小的值。

xgb_sk = XGBRegressor(max_depth=5,random_state=1412) #实例化
result_xgb_sk = cross_validate(xgb_sk,X,y,cv=cv
                               ,scoring="neg_root_mean_squared_error" #负根均方误差
                               ,return_train_score=True
                               ,verbose=True
                               ,n_jobs=-1)
def RMSE(result,name):
    return abs(result[name].mean())
train_score=RMSE(result_xgb_sk,"train_score")
test_score=RMSE(result_xgb_sk,"test_score")
print(f"train_score: {train_score:.2f}\ntest_score: {test_score:.2f}")
#train_score: 2992.28
#test_score: 27041.26

过拟合程度减轻，这说明模型是有潜力的，经过精密的调参之后xgboost上应该能够获得不错的结果。

xgb_sk = XGBRegressor(max_depth=5,random_state=12).fit(X,y)
xgb_sk.feature_importances_ #查看特征重要性

1.4.3 XGBoost回归原生代码实现

  原生代码必须使用XGBoost自定义的数据结构DMatrix，这一数据结构能够保证xgboost算法运行更快，并且能够自然迁移到GPU上运行。

import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import root_mean_squared_error

data_xgb = xgb.DMatrix(X,y) #1.将数据转换为DMatrix
x_train,x_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=12)
dtrain = xgb.DMatrix(x_train,y_train)
dtest = xgb.DMatrix(x_test,y_test)
params = {"max_depth":5,"seed":12} #2.定义参数
reg = xgb.train(params, data_xgb, num_boost_round=100) #3.训练模型
y_pred = reg.predict(data_xgb) #4.预测结果
root_mean_squared_error(y,y_pred) #RMSE #输出3296.846746

params = {"max_depth":5,"seed":12}
result = xgb.cv(params,data_xgb,num_boost_round=100
                ,nfold=5 #补充交叉验证中所需的参数，nfold=5表示5折交叉验证
                ,seed=12 #交叉验证的随机数种子，params中的是管理boosting过程的随机数种子
               )
result

import matplotlib.pyplot as plt
plt.figure(dpi=300)
plt.plot(result["train-rmse-mean"])
plt.plot(result["test-rmse-mean"])
plt.legend(["train","test"])
plt.title("xgboost 5fold cv")

1.5 XGBoost时间序列预测

  该数据为2006 年至 2016 年，匈牙利塞格德周边历史天气数据，一共96453条，包括：时间、天气类型、降水类型、气温、体表温度、湿度、风速（km/小时）、风向（度）、能见度（km）、噪声、气压、全天天气12个变量。

1.5.1 滑动窗口

  在时间序列预测问题中，滑动窗口是一种常用的数据处理方法，用于将时间序列数据转换为模型的输入特征和输出标签。滑动窗口的基本思想是以固定的时间窗口长度对时间序列进行切片，每次滑动一定的步长，从而生成一系列的子序列。这些子序列可以作为模型的输入特征，同时可以对应相同长度的下一个时间步的数据作为输出标签。这样就可以将时间序列数据转换为监督学习问题的数据集，用于训练和测试预测模型。
  通过滑动窗口的处理，原始的时间序列数据被转换为一系列的样本，每个样本包括了固定长度的输入特征和对应的输出标签，用于模型的训练和测试。滑动窗口技术可以帮助模型捕捉时间序列数据的局部模式和趋势，提高模型对时间序列的预测能力。

import numpy as np
 
# 定义参数
total_data_points = 16  # 总数据点数
window_size = 6         # 窗口大小
predict_length = 2      # 预测长度
 
# 随机生成一个数据集
data_set = np.random.rand(total_data_points)
# 分割数据集的函数
def create_inout_sequences(input_data, tw, pre_len):
    # 创建时间序列数据专用的数据分割器
    inout_seq = []
    L = len(input_data)
    for i in range(L - tw):
        train_seq = input_data[i:i + tw]
        if (i + tw + pre_len) > len(input_data):
            break
        train_label = input_data[i + tw:i + tw + pre_len]
        inout_seq.append((train_seq, train_label))
    return inout_seq
 
# 使用函数分割数据
X = create_inout_sequences(data_set, window_size, predict_length)
 
print(X)
#输出：[(('训练数据'),('标签')),(('训练数据'),('标签')),(('训练数据'),('标签'))

所得样本数=数据总数-(窗口大小+预测长度-1) --> 9=16-(6+2-1)
注意：在使用滑动窗口生成训练样本之前，通常需要对数据进行标准化或归一化，以提高模型的收敛速度和预测准确性。如果使用了归一化，输出预测结果时也要反归一化。

1.5.2 XGBoost 单变量预测

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.model_selection import TimeSeriesSplit
from sklearn.preprocessing import StandardScaler
from xgboost import XGBRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

def prepare_train_plot(df, lag_start=1, lag_end=5, test_size=0.2, target_encoding=False, plot_intervals=False, plot_anomalies=False):
    """
    1. 数据准备：创建滞后特征，时间特征，处理缺失值。
    2. 模型训练：使用XGBoost模型进行训练。
    3. 绘图：展示预测结果、实际结果，并可选地绘制置信区间与异常值。
    """
    # 数据准备
    data = pd.DataFrame(df.copy())
    data.columns = ["y"]  #将预测目标重命名为y，为y创建滞后特征
    for i in range(lag_start, lag_end + 1):
        data[f"lag_{i}"] = data.y.shift(i)
    
    # data["hour"] = data.index.hour  #小时
    data["weekday"] = data.index.weekday  #工作日
    data["is_weekend"] = (data["weekday"] >= 5).astype(int)  #周末
    #target_encoding参数为 True，则对 weekday 列进行目标编码（将不同星期几的目标变量均值作为新特征）
    if target_encoding:
        data["weekday_average"] = data["weekday"].map(data.groupby('weekday')['y'].mean())
        data.drop(["weekday"], axis=1, inplace=True) 
    
    data.dropna(inplace=True)
    y = data["y"]
    X = pd.get_dummies(data.drop(["y"], axis=1))
 
    # 数据分割
    test_index = int(len(X) * (1 - test_size))
    X_train, X_test = X.iloc[:test_index], X.iloc[test_index:]
    y_train, y_test = y.iloc[:test_index], y.iloc[test_index:]

    # 数据标准化
    scaler = StandardScaler()
    X_train_scaled = scaler.fit_transform(X_train)
    X_test_scaled = scaler.transform(X_test)

    # XGBoost训练
    xgb = XGBRegressor(n_estimators=100, learning_rate=0.05, random_state=42)
    xgb.fit(X_train_scaled, y_train)

    # 预测与绘图
    prediction = xgb.predict(X_test_scaled)
    plt.figure(figsize=(15, 7))
    plt.plot(prediction, "g", label="Prediction", linewidth=2.0)
    plt.plot(y_test.values, label="Actual", linewidth=2.0)

    # 置信区间
    if plot_intervals:
        tscv = TimeSeriesSplit(n_splits=5)
        cv = cross_val_score(xgb, X_train_scaled, y_train, cv=tscv, scoring="neg_mean_squared_error") #负均方误差
        
        '''cross_val_score 是 Scikit-learn 库中的一个用于评估模型性能的函数。它通过交叉验证来对模型进行训练和评估。
        scores = cross_val_score(model, X, y, cv=5, scoring='neg_mean_squared_error')
        cv：交叉验证的折数。默认为 5。也可以使用 StratifiedKFold 或 TimeSeriesSplit 来进行特定类型的交叉验证。
        scoring:用于评估模型性能的标准。分类如：accuracy, precision, recall, f1 等，回归：neg_mean_squared_error(负均方误差),r2(R²)'''
        
        deviation = np.sqrt(-cv.mean()) #均方根误差
        lower = prediction - (1.96 * deviation)
        upper = prediction + (1.96 * deviation)

        plt.fill_between(range(len(prediction)),lower, upper, color="r", alpha=0.2, label="Confidence Interval")

    # rmse = np.sqrt(mean_squared_error(y_test, prediction))
    r2 = r2_score(y_test, prediction)
    plt.title(f"R-squared: {r2:.2f}")
    plt.legend(loc="best")
    plt.grid(True)
    plt.tight_layout()
    plt.show()


if __name__ == '__main__':
     
    df = pd.read_csv('Temperature.csv', parse_dates=['date'])
    df.set_index('date', inplace=True)

    # 调用主函数：准备数据、训练模型并绘制结果
    prepare_train_plot(df[['Temperature']], lag_start=1, lag_end=5, test_size=0.2, target_encoding=True, plot_intervals=True, plot_anomalies=True)

plt.figure(figsize=(15,5), dpi=100)
plt.grid(True)
plt.plot(daily_avg_wind_speed['Temperature'], color='green')
plt.show()

1.5.3 XGBoost 多变量预测

from numpy import loadtxt
from xgboost import XGBClassifier, plot_tree
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt
import os
import numpy as np
import pandas as pd
import graphviz
from xgboost import to_graphviz

# 设置环境变量以避免一些错误
os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'

# 加载数据
dataset = loadtxt('pima-indians-diabetes.csv', delimiter=',')
X = dataset[:, 0:8]
y = dataset[:, 8]

# 将数据划分为训练集和测试集
seed = 7
test_size = 0.2
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=seed)

# 创建并训练模型
model = XGBClassifier(n_jobs=-1)
model.fit(X_train, y_train)

# 使用训练后的模型对测试集进行预测，并计算预测值与实际之间的准确率
y_pred = model.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print("Accuracy: %.2f%%" % (accuracy * 100.0))  #Accuracy: 74.03%

# 可视化XGBoost的树结构

# 选择树的索引，0表示第一棵树
tree_index = 0

# 使用plot_tree函数直接画出树
plot_tree(model, num_trees=tree_index)
plt.show()

# 或者使用to_graphviz将树转换为graphviz格式并渲染
# graph = to_graphviz(model, num_trees=tree_index)
# graph.render("xgboost_tree", format="png", cleanup=True)  # 将树保存为PNG文件

#使⽤训练后的模型对测试集进⾏预测，得到每个类别的预测概率
y_pred = model.predict(X_test)
y_pred_proba = model.predict_proba(X_test)
y_pred_proba
#输出各特征重要程度
from xgboost import plot_importance
from matplotlib import pyplot
%matplotlib inline
plot_importance(model)
pyplot.show()

#导⼊调参相关包
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import StratifiedKFold
#创建模型及参数搜索空间
model_GS = XGBClassifier()
learning_rate = [0.0001, 0.001, 0.01, 0.1, 0.2, 0.3]
max_depth = [1, 2, 3, 4, 5]
param_grid = dict(learning_rate=learning_rate, max_depth=max_depth)
#设置分层抽样验证及创建搜索对象
kfold = StratifiedKFold(n_splits=10, shuffle=True, random_state=seed)
grid_search = GridSearchCV(model_GS, param_grid=param_grid, scoring='neg_log_loss',n_jobs=-1, cv=kfold)
grid_result = grid_search.fit(X, y)      
y_pred = grid_result.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print("Accuracy: %.2f%%" % (accuracy * 100.0)) #Accuracy: 88.31%
grid_result.best_score_, grid_result.best_params_ #(np.float64(-0.4720604187971097), {'learning_rate': 0.1, 'max_depth': 2})

二、XGBoost特性总结

1.梯度提升：XGBoost采用梯度提升算法，通过逐步训练多个弱学习器（通常是决策树）来减少预测误差。每个新模型都是在前一模型的残差上进行优化，从而不断提升预测精度。
2.基于树的模型：XGBoost默认使用CART树作为基学习器，每棵树通过递归地划分特征空间来进行分类或回归。树模型能够捕捉数据中的复杂非线性关系，并且易于解释。
3.正则化策略：为防止过拟合，XGBoost引入了正则化手段，控制树的复杂度。通过限制树的最大深度、叶子节点的最小样本数和叶子节点的权重衰减等方式，有效控制模型的学习能力，增强其泛化性能。
4.特征选择和分裂：在构建决策树时，XGBoost通过选择最具信息增益的特征进行节点分裂，从而最小化预测误差。每次分裂都会寻找能带来最大增益的特征划分点，确保每一层的学习效果最优。
5. 并行计算：XGBoost通过并行计算显著提升训练速度。通过多线程和分布式计算，训练过程被分解成多个子任务并行执行，特别是在大规模数据集上表现出色，能够加速训练过程。
6. 自定义损失函数：XGBoost允许用户根据具体任务需求自定义损失函数，支持不同的回归或分类任务。通过灵活地定义损失函数，能够更好地适应各种实际应用场景。
Xgboost 防止过拟合的方法：

1. 目标函数添加正则化：叶子节点个数和叶子节点权重的正则化
2. 列抽样：训练的时候只用一部分特征
3. 子采样：每轮计算可以不使用全部样本，使算法更加保守