XGBoost实战详细代码——匈牙利水痘病例时间序列预测

目录

引言

数据集

XGBoost原理简述

代码及解析

结果展示

---

引言

        本文为机器学习大作业选题——匈牙利水痘病例预测（Prediction of the Hungarian case of chicken pox）的XGBoost解决方法，参与答辩并获得优秀，主要是使用XGBoost对病例数进行时间序列预测，包含数据集、XGBoost原理简述、代码及解析和结果展示四部分。所有文件如下

 下载链接：https://pan.baidu.com/s/1x3IAx_GFb2sBBznT59j7OA      提取码：1234 
 

数据集

        提供的数据集分为三部分，分别是包含422条时间序列数据的训练集、包含100条时间序列数据的测试集和包含各县关系的空间连接图，训练集和测试集记录了对应连续时间序列的匈牙利各县的每周患病人数，空间连接图记录了各县两两连接的关系，共有20个顶点（19个县和首都布达佩斯），节点之间有61条边。

XGBoost原理简述

        XGBoost本质上还是一个GBDT，但是力争把速度和效率发挥到极致，所以叫X (Extreme) GBoosted，GBDT(Gradient Boosting Decision Tree)，全名叫梯度提升决策树，使用的是Boosting的思想。Boosting方法训练基分类器时采用串行的方式，各个基分类器之间有依赖。它的基本思路是将基分类器层层叠加，每一层在训练的时候，对前一层基分类器分错的样本，给予更高的权重。测试时，根据各层分类器的结果的加权得到最终结果。

        基于Boosting，GBDT的原理就是所有弱分类器的结果相加等于预测值，然后下一个弱分类器去拟合误差函数对预测值的残差(这个残差就是预测值与真实值之间的误差)，它里面的弱分类器的表现形式就是CART。

        上述是我在阅读了一些博客和文章后的看法，下面是个人觉得对XGBoost和GDBT讲的比较好的博客链接

先对XGBoost有一个大致了解：终于有人说清楚了--XGBoost算法 - mantch

了解XGBoost绕不开GBDT：GBDT--原来是这么回事(附代码) - mantch

有了总体概念后可以深入了解细节：通俗理解kaggle比赛大杀器xgboost

这篇也很好，有示例和实战推荐：深入理解XGBoost，优缺点分析，原理推导及工程实现

        如果能看完并吃透上述四篇博客，基本对XGBoost就理解了，但是还是需要经过大量的实战操作才能真正掌握这门大杀器。

代码及解析

        代码的主体部分来自：如何用XGBoost做时间序列预测？ | 机器之心 (jiqizhixin.com)

        下文列出可运行代码及详细注释（添加了读取数据集、自动调优代码和可视化代码）

from numpy import asarray
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from sklearn.metrics import mean_absolute_error
from xgboost import XGBRegressor
from matplotlib import pyplot
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV

# transform a timeseries dataset into a supervised learning dataset
#这部分是对输入数据做转换，一列变成四列并错开，最好理解shift函数之后打印一下就懂了
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
    n_vars = 1 #if type(data) is list else data.shape[1]
    df = DataFrame(data)
    # print(df)
    cols = list()

    # input sequence (t-n, ... t-1)
    for i in range(n_in, 0, -1):
        cols.append(df.shift(i))

    # forecast sequence (t, t+1, ... t+n)
    for i in range(0, n_out):
        cols.append(df.shift(-i))

    # put it all together
    agg = concat(cols, axis=1)
    # print(agg)

    # drop rows with NaN values
    if dropnan:
        agg.dropna(inplace=True)
    return agg.values

# split a univariate dataset into train/test sets  分割训练集和测试集
def train_test_split(data, n_test):
    return data[:-n_test, :], data[-n_test:, :]


# fit an xgboost model and make a one step prediction
#训练模型并进行对测试集的预测，我加入了模型调优的代码
def xgboost_forecast(train, testX):
    # transform list into array
    train = asarray(train)

    # split into input and output columns
    trainX, trainy = train[:, :-1], train[:, -1]

    # 模型调优
    # cv_params = { 'n_estimators' : [200,250,300,350,400,450]}  #可以修改括号内的参数进行自动调参
    # other_params = { 'learning_rate': 0.1,  'n_estimators': 400,  'max_depth': 5,  'min_child_weight': 1,  'seed': 0,  'subsample': 0.8,  'colsample_bytree': 0.8,  'gamma': 0,  'reg_alpha': 0,  'reg_lambda': 1}
    # model = XGBRegressor(**other_params)
    # optimized_GBM = GridSearchCV(estimator=model, param_grid=cv_params, scoring='r2', cv=5, verbose=1, n_jobs=4)
    # optimized_GBM.fit(trainX, trainy)
    # evalute_result = optimized_GBM.cv_results_
    # print('每轮迭代运行结果:{0}' .format(evalute_result))
    # print('参数的最佳取值：{0}' .format(optimized_GBM.best_params_))
    # print('最佳模型得分:{0}' .format(optimized_GBM.best_score_))
    #训练模型
    model = XGBRegressor(objective='reg:squarederror', n_estimators=100)
    model.fit(trainX, trainy)

    # make a one-step prediction
    yhat = model.predict(asarray([testX]))
    return yhat[0]


# walk-forwardvalidation for univariate data
#封装分割数据集、每次添加一行数据进行训练和输出预测结果
def walk_forward_validation(data, n_test):
    predictions = list()

    # split dataset
    train, test = train_test_split(data, n_test)

    # seed history with training dataset
    history = [x for x in train]

    # step over each time-step in the test set
    for i in range(len(test)):
        # split test row into input and output columns
        testX, testy = test[i, :-1], test[i, -1]

        # fit model on history and make ap rediction
        yhat = xgboost_forecast(history, testX)

        # store forecast in list of predictions
        predictions.append(yhat)

        # add actual observation to history for the next loop
        history.append(test[i])

        # summarize progress
        print('>expected=%.1f,predicted=%.1f' % (testy, yhat))

    # estimate prediction error
    error = mean_absolute_error(test[:, 1], predictions)
    return error, test[:, 1], predictions

totalMSE = 0  #保存20个县的总体MSE

# load the dataset
series = read_csv('pre_522days_train_data.csv', header=0, index_col=0) #使用read_csv读取数据
values = series.values
# #对数据进行标准化处理，加了后拟合效果更好
# preprocess = StandardScaler()  
# values = preprocess.fit_transform(series)
# print(values)

for country in range(0,19):
    value = values[:, country]
    # pyplot.plot(values)
    # pyplot.show()  #可以先看一下原始数据

    # transform the timeseries data into supervised learning
    data = series_to_supervised(value, n_in=4)

    # evaluate
    mae, y, yhat = walk_forward_validation(data, 40)
    mse = 0
    for i in range(len(y)):
        mse = mse + (y[i] - yhat[i]) ** 2
        # print(mse)
    mse = mse / 100
    totalMSE = totalMSE + mse
    print('MAE: %.3f' % mae)
    print('MSE: %.3f' % mse)

    # plot expected vspreducted
    pyplot.plot(y, label='Expected')
    pyplot.plot(yhat, label='Predicted')
    pyplot.legend()
    pyplot.show()
totalMSE = totalMSE/20
print("totalMSE", totalMSE)

结果展示

        下图展示首都BUDAPEST和BARANYA、BACS、BEKES三县的预测结果，迭代次数、学习率和单棵CART最大深度分别设置为1000，0.2和5。

         数据标准化之后再进行预测，拟合效果会更好，MAE=0.508  MSE=0.444，下图为首都Budapest标准化之后的拟合效果，对比上图（a）效果好了一些