kaggle实战-精美可视化与时序预测

import math
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess
from sklearn.linear_model import LinearRegression
from pandas import date_range
from statsmodels.graphics.tsaplots import plot_pacf

from sklearn.model_selection import train_test_split

from sklearn.preprocessing import LabelEncoder
from xgboost import XGBRegressor

from sklearn.linear_model import ElasticNet, Lasso, Ridge

from sklearn.ensemble import ExtraTreesRegressor, RandomForestRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.multioutput import RegressorChain
import warnings
warnings.filterwarnings("ignore")

holidays = pd.read_csv('holidays_events.csv')
oil = pd.read_csv('oil.csv')
stores = pd.read_csv('stores.csv')
trans = pd.read_csv('transactions.csv')

train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')

# holidays

holidays.head()

holidays.shape

(350, 6)

holidays.describe()

holidays.isnull().sum()

date           0
type           0
locale         0
locale_name    0
description    0
transferred    0
dtype: int64

train.head()

oil.head()

stores.head()

trans.head()

holidays['date'] = pd.to_datetime(holidays['date'], format = "%Y-%m-%d")
oil['date'] = pd.to_datetime(oil['date'], format = "%Y-%m-%d")
trans['date'] = pd.to_datetime(trans['date'], format = "%Y-%m-%d")
train['date'] = pd.to_datetime(train['date'], format = "%Y-%m-%d")
test['date'] = pd.to_datetime(test['date'], format = "%Y-%m-%d")

oil.head()

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(25,15))
oil.plot.line(x="date",
              y="dcoilwtico",
              color='b',  # 颜色
              title ="dcoilwtico",  # 标题
              ax = axes,
              rot=0)
plt.show()

# 定义分组统计函数

def grouped(df, key, freq, col):
    """
    定义分组统计均值的函数
    """
    df_groupby = df.groupby([pd.Grouper(key=key,freq=freq)]).agg(mean=(col,'mean'))
    df_groupby = df_groupby.reset_index()
    return df_groupby

df_groupby_trans_w = grouped(trans, 'date', 'W', 'transactions')
df_groupby_trans_w.head(10)

def add_time(df, key, freq, col):
    # 调用上面的函数
    df_groupby = grouped(df, key, freq, col)
    # 添加time字段
    df_groupby['time'] = np.arange(len(df_groupby.index))
    # 将time字段放在索引为1的位置
    column_time = df_groupby.pop('time')
    df_groupby.insert(1, 'time', column_time)
    return df_groupby

# 基于week和month
df_groupby_train_w = add_time(train, 'date', 'W', 'sales')
df_groupby_train_m = add_time(train, 'date', 'M', 'sales')

df_groupby_train_w.head()

fig, axes = plt.subplots(nrows=3,
                         ncols=1,
                         figsize=(30,20)
                        )

# transactions基于weekly
axes[0].plot('date',
             'mean',
             data=df_groupby_trans_w,
             color='grey',
             marker='o')

axes[0].set_title("Transactions (groupby by week)", fontsize=20)

# sales-weekly
axes[1].plot('time',
             'mean',
             data=df_groupby_train_w,
             color='0.75')
axes[1].set_title("Sales (groupby by week)", fontsize=20)

# 线性回归
axes[1] = sns.regplot(x='time',
                      y='mean',
                      data=df_groupby_train_w,
                      scatter_kws=dict(color='0.75'),
                      ax = axes[1])

# sales-monthly
axes[2].plot('time',
             'mean',
             data=df_groupby_train_m,
             color='0.75')

axes[2].set_title("Sales (groupby by month)", fontsize=20)
# 基于线性回归
axes[2] = sns.regplot(x='time',
                      y='mean',
                      data=df_groupby_train_m,
                      scatter_kws=dict(color='0.75'),
                      line_kws={"color": "red"},
                      ax = axes[2])

plt.show()

def add_lag(df, key, freq, col, lag):
    """
    lag表示移动的标志位长度
    """
    df_groupby = grouped(df, key, freq, col)
    name = 'Lag_' + str(lag)
    df_groupby['Lag'] = df_groupby['mean'].shift(lag)
    return df_groupby

# 移动一个单位

df_groupby_train_w_lag1 = add_lag(train, 'date', 'W','sales',1)

df_groupby_train_m_lag1 = add_lag(train, 'date', 'W', 'sales', 1)

df_groupby_train_w_lag1.head()

fig, axes = plt.subplots(nrows=2,
                         ncols=1,
                         figsize=(30,20))

axes[0].plot('Lag',
             'mean',
             data=df_groupby_train_w_lag1,
             color='0.75',
             linestyle=(0, (1, 10)))

axes[0].set_title("Sales (groupby by week)", fontsize=20)
axes[0] = sns.regplot(x='Lag',
                      y='mean',
                      data=df_groupby_train_w_lag1,
                      scatter_kws=dict(color='0.75'),
                      ax = axes[0])


axes[1].plot('Lag',
             'mean',
             data=df_groupby_train_m_lag1,
             color='0.75',
             linestyle=(0, (1, 10)))

axes[1].set_title("Sales (groupby by month)", fontsize=20)
axes[1] = sns.regplot(x='Lag',
                      y='mean',
                      data=df_groupby_train_m_lag1,
                      scatter_kws=dict(color='0.75'),
                      line_kws={"color": "red"},
                      ax = axes[1])

plt.show()

def plot_stats(df, col, ax, color, angle):
    # 统计数量
    count_classes = df[col].value_counts()

    ax = sns.barplot(x=count_classes.index,
                     y=count_classes,
                     ax=ax,
                     palette=color)
    # 标题
    ax.set_title(col.upper(), fontsize=18)
    # 设置x轴旋转角度
    for tick in ax.get_xticklabels():
        tick.set_rotation(angle)

fig, axes = plt.subplots(nrows=1,
                         ncols=2,
                         figsize=(15,5))
fig.autofmt_xdate()
fig.suptitle("Stats of df_holidays".upper())
plot_stats(holidays, "type", axes[0], "pastel", 45)
plot_stats(holidays, "locale", axes[1], "rocket", 45)
plt.show()

fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(20,40))

plot_stats(stores, "city", axes[0], "mako_r", 45)
plot_stats(stores, "state", axes[1], "rocket_r", 45)
plot_stats(stores, "type", axes[2], "magma", 0)
plot_stats(stores, "cluster", axes[3], "viridis", 0)
plt.show()

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(20,10))
count_classes = train['family'].value_counts()
print("count_classes: \n", count_classes)

plt.title("Stats of train")
colors = ['#ff9999','#66b3ff','#99ff99',
          '#ffcc99', '#ffccf9', '#ff99f8',
          '#ff99af', '#ffe299', '#a8ff99',
          '#cc99ff', '#9e99ff', '#99c9ff',
          '#99f5ff', '#99ffe4', '#99ffaf']

plt.pie(count_classes,
        labels = count_classes.index,
        autopct='%1.1f%%',
        shadow=True,
        startangle=90,
        colors=colors)

plt.show()

def plot_boxplot(palette,x,y,hue,ax,title):
    # 设置主题
    sns.set_theme(style="ticks",palette=palette)
    # 绘图
    ax = sns.boxplot(x=x,y=y,hue=hue,ax=ax)
    # 设置标题和字体大小
    ax.set_title(title,fontsize=18)

fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10,20))

plot_boxplot("pastel", oil['date'].dt.year, oil['dcoilwtico'], oil['date'].dt.month, axes[0], "oil")
plot_boxplot("pastel", oil['date'].dt.year, oil['dcoilwtico'], oil['date'].dt.year, axes[1], "oil")
plt.show()

fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10,20))

plot_boxplot("pastel", trans['date'].dt.year, trans['transactions'], trans['date'].dt.month, axes[0], "trans")
plot_boxplot("pastel", trans['date'].dt.year, trans['transactions'], trans['date'].dt.year, axes[1], "trans")
plt.show()

# 0-定义分组统计函数-grouped

def grouped(df, key, freq, col):
    """
    定义分组统计均值的函数
    """
    df_groupby = df.groupby([pd.Grouper(key=key,freq=freq)]).agg(mean=(col,'mean'))
    df_groupby = df_groupby.reset_index()
    return df_groupby

# 1、groupby函数直接统计-以train数据为例

train.head()

# 2、直接统计：基于date和W进行统计sales的均值

train_groupby = train.groupby([pd.Grouper(key="date",freq="W")]).agg(mean=("sales",'mean'))
train_groupby = train_groupby.reset_index()
train_groupby

train_groupby.head(10)

# 3、加上移动平均功能
# center 表示以当前元素为中心进行统计
# min_periods表示窗口的最小观测值个数

# 一个窗口周期内至少要有4天的数据；前面3天不满足，均为NaN
train_groupby.rolling(window=7,min_periods=4).mean()

# 如何理解参数center

train_groupby.rolling(window=7,center=True,min_periods=4).mean()

# 3+中心+3

# 第一个元素前面没有数据，只能从本身向后取数据：0(中心) + 3----->[0:4]

train_groupby.iloc[:4,1].sum() / 4

192.27905005901843

# 第二个元素为中心：0+1(中心)+2/3/4----->[0:5]

train_groupby.iloc[0:5,1].sum() / 5

193.5360935220985

# 第三个元素为中心：0/1 + 2(中心)+3/4/5----->[0:6]

train_groupby.iloc[0:6, 1].sum() / 6

192.57810663810037

# 第四个元素为中心：0/1/2 + 3(中心)+4/5/6----->[0:7]

train_groupby.iloc[0:7, 1].sum() / 7

192.7887082607605

def plot_moving_average(df,key,freq,col,window,min_periods,ax,title):
    df_groupby = grouped(df,key,freq,col)
    moving_average = df_groupby["mean"].rolling(window=window,
                                                center=True,
                                                min_periods=min_periods).mean()
#     print("moving_average: \n", moving_average)

    ax = df_groupby["mean"].plot(color="0.75",linestyle="dashdot",ax=ax)
    ax = moving_average.plot(linewidth=3,color="g",ax=ax)
    ax.set_title(title,fontsize=18)

fig, axes = plt.subplots(nrows=2,
                         ncols=1,
                         figsize=(30,20))

# 基于train和trans数据
plot_moving_average(train, 'date', 'W', 'sales', 7, 4, axes[0], 'Sales Moving Average')
plot_moving_average(trans, 'date', 'W', 'transactions', 7, 4, axes[1], 'Transactions Moving Average')

plt.show()

def plot_deterministic_process(df, key, freq, col, ax1, title1, ax2, title2):
    # 分组统计
    df_groupby = grouped(df, key, freq, col)
    df_groupby["date"] = pd.to_datetime(df_groupby["date"],format="%Y-%m-%d")

    dp = DeterministicProcess(index=df_groupby["date"],
                              constant=True,
                              order=1,
                              drop=True)
    # 将频率设置成索引
    dp.index.freq = freq

    X1 = dp.in_sample()   # 属性特征X
    y1 = df_groupby["mean"]  # 待预测的目标值y
    y1.index = X1.index
    # 实例化类
    model = LinearRegression(fit_intercept=False)
    model.fit(X1, y1)

    # 预测结果
    y1_pred = pd.Series(model.predict(X1), index=X1.index)

    ax1 = y1.plot(linestyle='dashed', label="mean", color="0.75", ax=ax1, use_index=True)
    ax1 = y1_pred.plot(linewidth=3, label="Trend", color='r', ax=ax1, use_index=True)
    # 设置标题
    ax1.set_title(title1, fontsize=18)
    _ = ax1.legend()


    # 预测未来30个步长的数据
    steps = 30
    X2 = dp.out_of_sample(steps=steps)
    y2_fore = pd.Series(model.predict(X2), index=X2.index)

    ax2 = y1.plot(linestyle='dashed', label="mean", color="0.75", ax=ax2, use_index=True)
    ax2 = y1_pred.plot(linewidth=3, label="Trend", color='r', ax=ax2, use_index=True)
    ax2 = y2_fore.plot(linewidth=3, label="Predicted Trend", color='r', ax=ax2, use_index=True)
    ax2.set_title(title2, fontsize=18)
    _ = ax2.legend()

fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(30,30))

plot_deterministic_process(train, 'date', 'W', 'sales',
                           axes[0], "Sales Linear Trend",
                           axes[1], "Sales Linear Trend Forecast")
plot_deterministic_process(trans, 'date', 'W', 'transactions',
                           axes[2], "Transactions Linear Trend",
                           axes[3], "Transactions Linear Trend Forecast")

plt.show()

def seasonal_plot(X,y,period,freq,ax=None):
    """
    seasonal_plot函数
    """
    if ax is None:
        _, ax = plot.subplots()

    # 颜色设置
    palette = sns.color_palette("husl", n_colors=X[period].nunique(),)

    ax = sns.lineplot(x=X[freq],
                     y=X[y],
                     ax=ax,
                     hue=X[period],
                     palette=palette,
                     legend=False)
    # 标题设置
    ax.set_title(f'Seasonal Plot({period}/{freq})')

    for line, name in zip(ax.lines, X[period].unique()):
        y_ = line.get_ydata()[-1]
        ax.annotate(name,
                   xy=(1,y_),
                   xytext=(6,0),
                   color=line.get_color(),
                   xycoords=ax.get_yaxis_transform(),
                   textcoords="offset points",
                   size=14,
                   va="center")

    return ax

def plot_periodogram(ts, detrend="linear", ax=None):
    """
    plot_periodogram函数：基于周期的直方图
    """

    from scipy.signal import periodogram  # 导入模块
    fs = pd.Timedelta("365D") / pd.Timedelta("1D")

    freqencies, spectrum = periodogram(ts,
                                      fs=fs,
                                      detrend=detrend,
                                      window="boxcar",
                                      scaling="spectrum")

    if ax is None:
        _, ax = plt.subplots()

    ax.step(freqencies, spectrum, color="purple")
    ax.set_xscale("log")
    ax.set_xticks([1, 2, 4, 6, 12, 26, 52, 104])
    ax.set_xticklabels(["Annual (1)", "Semiannual (2)", "Quarterly (4)",
                        "Bimonthly (6)", "Monthly (12)", "Biweekly (26)",
                        "Weekly (52)", "Semiweekly (104)"], rotation=35)

    ax.ticklabel_format(axis="y", style="sci", scilimits=(0, 0))
    ax.set_ylabel("Variance")
    ax.set_title("Periodogram")
    return ax

def seasonality(df, key, freq, col):
    """
    seasonality函数：调用seasonal_plot 和 plot_periodogram函数
    """
    df_groupby = grouped(df, key, freq, col) # 调用grouped函数
    df_groupby['date'] = pd.to_datetime(df_groupby['date'], format = "%Y-%m-%d")
    df_groupby.index = df_groupby['date']  # 索引设置
    df_groupby = df_groupby.drop(columns=['date'])  # 删除date列
    df_groupby.index.freq = freq #  索引设置成频率

    X = df_groupby.copy()  # 副本
    X.index = pd.to_datetime(X.index, format = "%Y-%m-%d")
    X.index.freq = freq  # 索引设置成频率

    # 一周中的第几天
    X["day"] = X.index.dayofweek
    X["week"] = pd.Int64Index(X.index.isocalendar().week)
    # 一年中的第几天
    X["dayofyear"] = X.index.dayofyear
    X["year"] = X.index.year

    fig, (ax0, ax1, ax2) = plt.subplots(3, 1, figsize=(20, 30))

    # 调用seasonal_plot函数
    seasonal_plot(X, y='mean', period="week", freq="day", ax=ax0)  # 周期按照week
    seasonal_plot(X, y='mean', period="year", freq="dayofyear", ax=ax1)  # 周期按照year
    X_new = (X['mean'].copy()).dropna()

    # 调用plot_periodogram函数
    plot_periodogram(X_new, ax=ax2)

seasonality(trans, 'date', 'D', 'transactions')

seasonality(train, 'date', 'D', 'sales')

def predict_seasonality(df, key, freq, col, ax1, title1):
    fourier = CalendarFourier(freq="A", order=10)
    df_groupby = grouped(df,key,freq,col)

    # 时间格式转化
    df_groupby['date'] = pd.to_datetime(df_groupby['date'], format = "%Y-%m-%d")
    # 频率设置
    df_groupby['date'].freq = freq

    dp = DeterministicProcess(index=df_groupby["date"],
                             constant=True,
                             order=1,
                             period=None,
                             seasonal=True,
                             additional_terms=[fourier],
                             drop=True)

    dp.index.freq = freq

    # X-y
    X1 = dp.in_sample()
    y1 = df_groupby["mean"]
    y1.index = X1.index

    # 线性回归模型实例化
    model = LinearRegression(fit_intercept=False)
    model.fit(X1, y1)
    # 模型预测
    y1_pred = pd.Series(model.predict(X1), index=X1.index)
    X1_fore = dp.out_of_sample(steps=90)
    y1_fore = pd.Series(model.predict(X1_fore), index=X1_fore.index)
    # 绘图
    ax1 = y1.plot(linestyle='dashed', style='.', label="init mean values", color="0.4", ax=ax1, use_index=True)
    ax1 = y1_pred.plot(linewidth=3, label="Seasonal", color='b', ax=ax1, use_index=True)
    ax1 = y1_fore.plot(linewidth=3, label="Seasonal Forecast", color='r', ax=ax1, use_index=True)
    # 标题
    ax1.set_title(title1, fontsize=18)
    _ = ax1.legend()

fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(30,30))

predict_seasonality(trans, 'date', 'W', 'transactions', axes[0], "Transactions Seasonal Forecast")
predict_seasonality(train, 'date', 'W', 'sales', axes[1], "Sales Seasonal Forecast")
plt.show()

store_sales = train.copy() # copy副本

store_sales["date"] = store_sales.date.dt.to_period("D")  # 提取日期D
store_sales = store_sales.set_index(["store_nbr","family","date"]).sort_index()  # 设置双层索引
store_sales.head()

family_sales = (
    store_sales  # 数据df
    .groupby(['family', 'date'])  # 分组对象
    .mean()   # 统计函数
    .unstack('family')  # 反堆叠对象
    .loc['2017', ['sales', 'onpromotion']]  # 位置索引
)

family_sales.head()

mag_sales = family_sales.loc(axis=1)[:, 'MAGAZINES']
mag_sales.head(10)

y = mag_sales.loc[:,"sales"].squeeze()
y

date
2017-01-01    0.074074
2017-01-02    7.777778
2017-01-03    3.500000
2017-01-04    3.500000
2017-01-05    3.203704
                ...
2017-08-11    9.259259
2017-08-12    8.944444
2017-08-13    8.685185
2017-08-14    8.462963
2017-08-15    8.537037
Freq: D, Name: MAGAZINES, Length: 227, dtype: float64

fourier = CalendarFourier(freq="M", order=4)
dp = DeterministicProcess(constant=True,
                          index=y.index,
                          order=1,
                          seasonal=True,
                          drop=True,
                          additional_terms=[fourier]
                         )

X_time = dp.in_sample()
# 元旦
X_time["NewYearsDay"] = (X_time.index.dayofyear == 1)

model = LinearRegression(fit_intercept=False)
model.fit(X_time,y)

y_deseason = y - model.predict(X_time)
y_deseason.name = "sales_deseasoned"

ax = y_deseason.plot()
ax.set_title("Magazine Sales (deseasonalized)")

Text(0.5, 1.0, 'Magazine Sales (deseasonalized)')

def lagplot(x, y=None, lag=1, standardize=False, ax=None, **kwargs):
    """
    lagplot函数
    """
    from matplotlib.offsetbox import AnchoredText

    x_ = x.shift(lag)  # 移位函数
    if standardize:  # 数据标准化
        x_ = (x_ - x_.mean()) / x_.std()
    if y is not None:
        y_ = (y - y.mean()) / y.std() if standardize else y
    else:
        y_ = x

    # 相关系数
    corr = y_.corr(x_)
    if ax is None:
        fig, ax = plt.subplots()
    scatter_kws = dict(alpha=0.75,s=3)
    line_kws = dict(color='C3', )
    ax = sns.regplot(x=x_,
                     y=y_,
                     scatter_kws=scatter_kws,
                     line_kws=line_kws,
                     lowess=True,
                     ax=ax,
                     **kwargs)

    at = AnchoredText(f"{corr:.2f}",
                      prop=dict(size="large"),
                      frameon=True,
                      loc="upper left"
                     )

    at.patch.set_boxstyle("square, pad=0.0")
    ax.add_artist(at)
    ax.set(title=f"Lag {lag}", xlabel=x_.name, ylabel=y_.name)
    return ax

def plot_lags(x, y=None, lags=6, nrows=1, lagplot_kwargs={}, **kwargs):
    import math
    kwargs.setdefault('nrows', nrows)
    kwargs.setdefault('ncols', math.ceil(lags / nrows))
    kwargs.setdefault('figsize', (kwargs['ncols'] * 2 + 10, nrows * 2 + 5))
    fig, axs = plt.subplots(sharex=True,
                            sharey=True,
                            squeeze=False,
                            **kwargs)

    for ax, k in zip(fig.get_axes(), range(kwargs['nrows'] * kwargs['ncols'])):
        if k + 1 <= lags:
            ax = lagplot(x,
                         y,
                         lag=k + 1,
                         ax=ax,
                         **lagplot_kwargs)

            ax.set_title(f"Lag {k + 1}",
                         fontdict=dict(fontsize=14))
            ax.set(xlabel="", ylabel="")
        else:
            ax.axis('off')
    plt.setp(axs[-1, :], xlabel=x.name)
    plt.setp(axs[:, 0], ylabel=y.name if y is not None else x.name)
    fig.tight_layout(w_pad=0.1, h_pad=0.1)
    return fig

fig = plot_lags(y_deseason, lags=8, nrows=2)  # 多行多列数据

fig = plot_pacf(y_deseason, lags=8)  # 一个图形

onpromotion = mag_sales.loc[:,'onpromotion'].squeeze().rename('onpromotion')

# 去除离群点：每年的第一天
plot_lags(x=onpromotion.iloc[1:], y=y_deseason.iloc[1:], lags=3, nrows=1)

def make_lags(ts, lags):
    return pd.concat(
        {
            f'y_lag_{i}': ts.shift(i)
            for i in range(1, lags + 1)
        },
        axis=1)

X = make_lags(y_deseason, lags=4)
X = X.fillna(0.0)

X

y = y_deseason.copy()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=60, shuffle=False)

# 线性回归
model = LinearRegression()
# 模型训练
model.fit(X_train, y_train)
# 预测
y_pred = pd.Series(model.predict(X_train), index=y_train.index)
y_fore = pd.Series(model.predict(X_test), index=y_test.index)

# 绘图
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(20,10))
ax = y_train.plot(color="0.75",
                  style=".-",
                  markeredgecolor="0.25",
                  markerfacecolor="0.25",
                  ax=ax)

ax = y_test.plot(color="0.75",
                 style=".-",
                 markeredgecolor="0.25",
                 markerfacecolor="0.25",
                 ax=ax)

ax = y_pred.plot(ax=ax)
_ = y_fore.plot(ax=ax, color='C3')

plt.show()

store_sales = train.copy()

store_sales['date'] = store_sales.date.dt.to_period('D')

store_sales = store_sales.set_index(['store_nbr', 'family', 'date']).sort_index()

family_sales = (
    store_sales
    .groupby(['family', 'date'])
    .mean()
    .unstack('family')
    .loc['2017']
)
family_sales.head()

class BoostedHybrid:
    def __init__(self, model_1, model_2):
        self.model_1 = model_1
        self.model_2 = model_2
        self.y_columns = None


def fit(self, X_1, X_2, y):
    """
    模型训练fit
    """
    # 训练模型1
    self.model_1.fit(X_1, y)

    # 预测
    y_fit = pd.DataFrame(
        self.model_1.predict(X_1),
        index=X_1.index,
        columns=y.columns
    )

    # 计算残差
    y_resid = y - y_fit
    # 宽表转成长表
    y_resid = y_resid.stack().squeeze()

    # 基于残差训练模型2
    self.model_2.fit(X_2, y_resid)
    self.y_columns = y.columns
    self.y_fit = y_fit
    self.y_resid = y_resid

# 将fit方法添加类中
BoostedHybrid.fit = fit

def predict(self, X_1, X_2):
    """
    模型预测
    """
    y_pred = pd.DataFrame(
        self.model_1.predict(X_1),
        index=X_1.index,
        columns=self.y_columns
    )

    # 宽表转成长表
    y_pred = y_pred.stack().squeeze()
    # 将两个模型预测的结果累加
    y_pred += self.model_2.predict(X_2)

    return y_pred.unstack()

# 将预测predict添加到类中
BoostedHybrid.predict = predict

# 目标值y
y = family_sales.loc[:, 'sales']

dp = DeterministicProcess(index=y.index, order=1)
X_1 = dp.in_sample()

# X_2：XBoost
X_2 = family_sales.drop('sales', axis=1).stack()

# family类比编码
le = LabelEncoder()
X_2 = X_2.reset_index('family')
X_2['family'] = le.fit_transform(X_2['family'])

# 季节性特征
X_2["day"] = X_2.index.day

model = BoostedHybrid(
    model_1 = LinearRegression(),
    model_2 = XGBRegressor()
)

model.fit(X_1, X_2, y)

y_pred = model.predict(X_1,X_2)
y_pred = y_pred.clip(0.0)

# 模型混合

model = BoostedHybrid(
    model_1=Ridge(),
    model_2=KNeighborsRegressor(),
)

# 构造训练集和验证集

y_train, y_valid = y[:"2017-07-01"], y["2017-07-02":]
X1_train, X1_valid = X_1[: "2017-07-01"], X_1["2017-07-02" :]
X2_train, X2_valid = X_2.loc[:"2017-07-01"], X_2.loc["2017-07-02":]

# 训练
model.fit(X1_train, X2_train, y_train)

# 预测
y_fit = model.predict(X1_train, X2_train).clip(0.0)
y_pred = model.predict(X1_valid, X2_valid).clip(0.0)

families = y.columns[0:6]
axs = y.loc(axis=1)[families].plot(subplots=True,
                                   sharex=True,
                                   figsize=(30, 20),
                                   color="0.75",
                                   style=".-",
                                   markeredgecolor="0.25",
                                   markerfacecolor="0.25",
                                   alpha=0.5)
_ = y_fit.loc(axis=1)[families].plot(subplots=True,
                                     sharex=True,
                                     color='C0',
                                     ax=axs)
_ = y_pred.loc(axis=1)[families].plot(subplots=True,
                                      sharex=True,
                                      color='C3',
                                      ax=axs)
for ax, family in zip(axs, families):
    ax.legend([])
    ax.set_ylabel(family)

store_sales = train.copy()
store_sales['date'] = store_sales.date.dt.to_period('D')
store_sales = store_sales.set_index(['store_nbr', 'family', 'date']).sort_index()

family_sales = (
    store_sales
    .groupby(['family', 'date'])
    .mean()
    .unstack('family')
    .loc['2017']
)

test = test.copy()
test['date'] = test.date.dt.to_period('D')
test = test.set_index(['store_nbr', 'family', 'date']).sort_index()

def make_multistep_target(ts, steps):
    return pd.concat(
        {f'y_step_{i + 1}': ts.shift(-i)
         for i in range(steps)},
        axis=1)

y = family_sales.loc[:, 'sales']

X = make_lags(y, lags=5).dropna()

# 构建多步预测
y = make_multistep_target(y, steps=16).dropna()

y, X = y.align(X, join='inner', axis=0)

# 准备好输入到XGBoost模型的数据

le = LabelEncoder()
X = (X
    .stack('family')
    .reset_index('family')
    .assign(family=lambda x: le.fit_transform(x.family))  # 编码
)
y = y.stack('family')  # 宽表转成长表

y

# 模型实例化
model = RegressorChain(base_estimator=XGBRegressor())

# 模型预测
model.fit(X, y)
# 预测结果
y_pred = pd.DataFrame(model.predict(X),
                      index=y.index,
                      columns=y.columns).clip(0.0)

def plot_multistep(y, every=1, ax=None, palette_kwargs=None):
    palette_kwargs_ = dict(palette='husl', n_colors=16, desat=None)
    if palette_kwargs is not None:
        palette_kwargs_.update(palette_kwargs)
    palette = sns.color_palette(**palette_kwargs_)
    if ax is None:
        fig, ax = plt.subplots()
    ax.set_prop_cycle(plt.cycler('color', palette))
    for date, preds in y[::every].iterrows():
        preds.index = pd.period_range(start=date, periods=len(preds))
        preds.plot(ax=ax)
    return ax

FAMILY = 'BEAUTY'
START = '2017-04-01'
EVERY = 16

y_pred_ = y_pred.xs(FAMILY, level='family', axis=0).loc[START:]
y_ = family_sales.loc[START:, 'sales'].loc[:, FAMILY]

fig, ax = plt.subplots(1, 1, figsize=(11, 4))
ax = y_.plot(color="0.75",style=".-",markeredgecolor="0.25", markerfacecolor="0.25",ax=ax, alpha=0.5)
ax = plot_multistep(y_pred_, ax=ax, every=EVERY)
_ = ax.legend([FAMILY, FAMILY + ' Forecast'])