Kaggle M5 Forecasting：傳統預測方法與機器學習預測方法對比

本文使用的數據集來自 kaggle：M5 Forecasting — Accuracy。 該數據集包含有 California、Texas、Wisconsin 三個州的產品類別、部門、倉儲信息等。基於這些數據，需要預測接下來 28 天的每日銷售量。

涉及到的方法有：

· 單指數平滑法

· 雙指數平滑法

· 三指數平滑法

· ARIMA

· SARIMA

· SARIMAX

· Light Gradient Boosting

· Random Forest

· Linear Regression

為了使用上述方法，首先導入相應的包/庫：

<code>import time
import warnings
import numpy as np
import pandas as pd 
import seaborn as sns 
import lightgbm as lgb
from itertools import cycle
from sklearn.svm import SVR 

import statsmodels.api as sm
from pmdarima import auto_arima
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.tsa.holtwinters import SimpleExpSmoothing, ExponentialSmoothing

%matplotlib inline 
plt.style.use('bmh')
sns.set_style("whitegrid")
plt.rc('xtick', labelsize=15) 
plt.rc('ytick', labelsize=15) 
warnings.filterwarnings("ignore")
pd.set_option('max_colwidth', 100)
pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)
color_pal = plt.rcParams['axes.prop_cycle'].by_key()['color']
color_cycle = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])/<code>

然後導入數據集：

<code>data = pd.read_csv('data_for_tsa.csv')data['date'] = pd.to_datetime(data['date'])data.head()/<code>

數據集前五行數據

數據集包含了 2011-01-29 到 2016-05-22 期間的 1941 天的數據。其中最後 28 天作為測試集。

預測目標是 demand，即：當日的產品銷售量。

接下來進行數據集劃分

測試集包含了 2016-03-27 到 2016-04-24 期間的 28 天的數據。 2016-03-27 之前的其他數據則作為訓練數據。

<code>train = data[data['date'] <= '2016-03-27']
test = data[(data['date'] > '2016-03-27') & (data['date'] <= '2016-04-24')]

fig, ax = plt.subplots(figsize=(25,5))
train.plot(x='date',y='demand',label='Train',ax=ax)
test.plot(x='date',y='demand',label='Test',ax=ax);/<code>

時間序列數據

為了便於對比所有方法的準確性，建立一個命名為 predictions 的 dataframe，將每個方法設為其中的一行。建立一個命名為stats的 dataframe，用於存儲每個方法的性能表現和計算時間。

<code>predictions = pd.DataFrame()
predictions['date'] = test['date']
stats = pd.DataFrame(columns=['Model Name','Execution Time','RMSE'])/<code>

訓練及評價模型

單指數平滑方法

通過調用 SimpleExpSmoothing 包，可以使用 EWMA, Exponentially Weighted Moving Average方法——一種單指數平滑方法。

使用 EWMA 方法，我們首先需要定義 span 變量——數據集的季節週期。

<code>fig, ax = plt.subplots(figsize=(15, 3))
plot_acf(data['demand'].tolist(), lags=60, ax=ax);/<code>

自相關圖

查看數據的自相關圖可知，每隔七個數據，達到一個峰值，也就意味著任一數據與之前的第七個時間數據具有較高的相關性。所以這裡將 span 設為 7。

具體地，通過以下代碼實現單指數平滑方法預測：

<code>t0 = time.time()
model_name='Simple Exponential Smoothing'
span = 7
alpha = 2/(span+1)
#train
simpleExpSmooth_model = SimpleExpSmoothing(train['demand']).fit(smoothing_level=alpha,optimized=False)
t1 = time.time()-t0
#predict
predictions[model_name] = simpleExpSmooth_model.forecast(28).values
#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t1, 'RMSE':score},ignore_index=True)/<code>

單指數平滑方法預測結果

上述代碼實現了對於數據的學習，通過 forcast(x),x=28，預測了接下來 28 天的數據。並且通過均方根誤差衡量誤差。

雙指數平滑方法

單指數平滑方法只使用了一個平滑係數 ，而雙指數平滑方法則引入了第二個平滑係數 ，以反映數據的趨勢。 使用雙指數平滑方法，我們需要定義 seasonal_periods。

具體代碼如下：

<code>t0 = time.time()
model_name='Double Exponential Smoothing'
#train
doubleExpSmooth_model = ExponentialSmoothing(train['demand'],trend='add',seasonal_periods=7).fit()
t1 = time.time()-t0
#predict
predictions[model_name] = doubleExpSmooth_model.forecast(28).values
#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t1, 'RMSE':score},ignore_index=True)/<code>

雙指數平滑方法預測結果

三指數平滑方法

三指數平滑方法進一步引入了係數以反映數據的趨勢及季節性變化。

具體代碼如下：

<code>t0 = time.time()
model_name='Triple Exponential Smoothing'
#train
tripleExpSmooth_model = ExponentialSmoothing(train['demand'],trend='add',seasonal='add',seasonal_periods=7).fit()
t1 = time.time()-t0
#predict
predictions[model_name] = tripleExpSmooth_model.forecast(28).values
#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t1, 'RMSE':score},ignore_index=True)/<code>

三指數平滑方法預測結果

從預測結果可以看出，三指數平滑方法能夠學習數據的季節性變化特徵。

ARIMA

使用 ARIMA 方法，首先需要確定 p,d,q 三個參數。

· p 是AR項的順序。

· d 是使時間序列平穩所需的差分次數

· q 是MA項的順序。

自動確定 ARIMA 所需參數

通過調用 auto_arima 包，可以自動確定 ARIMA 所需的參數。

t0 = time.time()model_name='ARIMA'arima_model = auto_arima(train['demand'], start_p=0, start_q=0, max_p=14, max_q=3, seasonal=False, d=None, trace=True,random_state=2020, error_action='ignore', # we don't want to know if an order does not work suppress_warnings=True, # we don't want convergence warnings stepwise=True)arima_model.summary()

auto_arima 的計算結果

確定了 p,d,q 參數，就可以進行下一步的訓練及預測：

<code>#train
arima_model.fit(train['demand'])
t1 = time.time()-t0
#predict
predictions[model_name] = arima_model.predict(n_periods=28)
#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t1, 'RMSE':score},ignore_index=True)/<code>

ARIMA 預測結果

這裡使用的簡單ARIMA模型不考慮季節性，是一個（5，1，3）模型。這意味著它使用5個滯後來預測當前值。移動窗口的大小等於 1，即滯後預測誤差的數量等於1。使時間序列平穩所需的差分次數為 3。

SARIMA

SARIMA 是 ARIMA 的發展，進一步引入了相關參數以使得模型能夠反映數據的季節變化特徵。

通過 auto_arima 相關代碼，將參數設置為 seasonal=True,m=7，自動計算 SARIMA 所需的參數。

<code>t0 = time.time()
model_name='SARIMA'
sarima_model = auto_arima(train['demand'], start_p=0, start_q=0,
                          max_p=14, max_q=3,
                          seasonal=True, m=7,
                          d=None, trace=True,random_state=2020,
                          out_of_sample_size=28,
                          error_action='ignore',   # we don't want to know if an order does not work
                          suppress_warnings=True,  # we don't want convergence warnings
                          stepwise=True)
sarima_model.summary() /<code>

auto_arima 計算結果

確定了參數後，接下來進行訓練及預測：

<code>#train
sarima_model.fit(train['demand'])
t1 = time.time()-t0
#predict
predictions[model_name] = sarima_model.predict(n_periods=28)
#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t1, 'RMSE':score},ignore_index=True)/<code>

SARIMA 預測結果

SARIMAX

使用前面的方法，我們只能基於前面的歷史數據進行預測。在 SARIMAX 中引入外生迴歸因子（eXogenous regressors），可以實現對時間序列數據以外的數據的分析。

本例中，我們引入 sell_price 數據以輔助更好地預測。

<code>t0 = time.time()
model_name='SARIMAX'
sarimax_model = auto_arima(train['demand'], start_p=0, start_q=0,
                          max_p=14, max_q=3,
                          seasonal=True, m=7,
                          exogenous = train[['sell_price']].values,
                          d=None, trace=True,random_state=2020,
                          out_of_sample_size=28,
                          error_action='ignore',   # we don't want to know if an order does not work
                          suppress_warnings=True,  # we don't want convergence warnings
                          stepwise=True)
sarimax_model.summary() /<code>

auto_arima 計算結果

通過 auto_arima 自動計算了 SARIMAX 方法所需的參數後，可以直接進行訓練和預測。

<code>#train
sarimax_model.fit(train['demand'])
t1 = time.time()-t0
#predict
predictions[model_name] = sarimax_model.predict(n_periods=28)
#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t1, 'RMSE':score},ignore_index=True)/<code>

SARIMA 預測結果

從預測結果可以看出，通過分析額外的數據，有助於減少誤差。

機器學習

使用機器學習方法，首先需要特徵數據以及指標數據。

在本文中，基於時間序列數據構造特徵數據如下：

· 特徵數據1：滯後數據。 選擇 7 天前的 demand 數據作為特徵數據。

· 特徵數據2：移動平均數據。選擇 7 天前至 14 天之前的 demand 移動平均值數據作為特徵數據。

· 特徵數據3：月銷售均值

· 特徵數據4：每月銷售最大值

· 特徵數據5：每月銷售最小值

· 特徵數據6：每月銷售最大值與最小值的差值

· 特徵數據7：每週銷售均值

· 特徵數據8：每週銷售最大值

· 特徵數據9：每週銷售中值

具體代碼如下：

<code>def lags_windows(df):
    lags = [7]
    lag_cols = ["lag_{}".format(lag) for lag in lags ]
    for lag, lag_col in zip(lags, lag_cols):
        df[lag_col] = df[["id","demand"]].groupby("id")["demand"].shift(lag)

    wins = [7]
    for win in wins :
        for lag,lag_col in zip(lags, lag_cols):
            df["rmean_{}_{}".format(lag,win)] = df[["id", lag_col]].groupby("id")[lag_col].transform(lambda x : x.rolling(win).mean())   
    return df

def per_timeframe_stats(df, col):
    #For each item compute its mean and other descriptive statistics for each month and dayofweek in the dataset
    months = df['month'].unique().tolist()
    for y in months:
        df.loc[df['month'] == y, col+'_month_mean'] = df.loc[df['month'] == y].groupby(['id'])[col].transform(lambda x: x.mean()).astype("float32")
        df.loc[df['month'] == y, col+'_month_max'] = df.loc[df['month'] == y].groupby(['id'])[col].transform(lambda x: x.max()).astype("float32")
        df.loc[df['month'] == y, col+'_month_min'] = df.loc[df['month'] == y].groupby(['id'])[col].transform(lambda x: x.min()).astype("float32")
        df[col + 'month_max_to_min_diff'] = (df[col + '_month_max'] - df[col + '_month_min']).astype("float32")

    dayofweek = df['dayofweek'].unique().tolist()
    for y in dayofweek:
        df.loc[df['dayofweek'] == y, col+'_dayofweek_mean'] = df.loc[df['dayofweek'] == y].groupby(['id'])[col].transform(lambda x: x.mean()).astype("float32")
        df.loc[df['dayofweek'] == y, col+'_dayofweek_median'] = df.loc[df['dayofweek'] == y].groupby(['id'])[col].transform(lambda x: x.median()).astype("float32")
        df.loc[df['dayofweek'] == y, col+'_dayofweek_max'] = df.loc[df['dayofweek'] == y].groupby(['id'])[col].transform(lambda x: x.max()).astype("float32")
    return df

def feat_eng(df):
    df = lags_windows(df)
    df = per_timeframe_stats(df,'demand')
    return df/<code>

準備數據：

<code>data = pd.read_csv('data_for_tsa.csv')
data['date'] = pd.to_datetime(data['date'])
train = data[data['date'] <= '2016-03-27']
test = data[(data['date'] > '2016-03-11') & (data['date'] <= '2016-04-24')]
  
data_ml = feat_eng(train)
data_ml = data_ml.dropna() 

useless_cols = ['id','item_id','dept_id','cat_id','store_id','state_id','demand','date','demand_month_min']
linreg_train_cols = ['sell_price','year','month','dayofweek','lag_7','rmean_7_7'] #use different columns for linear regression
lgb_train_cols = data_ml.columns[~data_ml.columns.isin(useless_cols)] 

X_train = data_ml[lgb_train_cols].copy()
y_train = data_ml["demand"]/<code>

模型擬合

通過 light gradient boosting、linear regression、random forest 三種方法對數據進行擬合：

<code>#Fit Light Gradient Boosting
t0 = time.time()
lgb_params = {
        "objective" : "poisson",
        "metric" :"rmse",
        "force_row_wise" : True,
        "learning_rate" : 0.075,
        "sub_row" : 0.75,
        "bagging_freq" : 1,
        "lambda_l2" : 0.1,
        'verbosity': 1,
        'num_iterations' : 2000,
        'num_leaves': 128,
        "min_data_in_leaf": 50,
}
np.random.seed(777)
fake_valid_inds = np.random.choice(X_train.index.values, 365, replace = False)
train_inds = np.setdiff1d(X_train.index.values, fake_valid_inds)
train_data = lgb.Dataset(X_train.loc[train_inds] , label = y_train.loc[train_inds], free_raw_data=False)
fake_valid_data = lgb.Dataset(X_train.loc[fake_valid_inds], label = y_train.loc[fake_valid_inds],free_raw_data=False)

m_lgb = lgb.train(lgb_params, train_data, valid_sets = [fake_valid_data], verbose_eval=0) 
t_lgb = time.time()-t0

#Fit Linear Regression
t0 = time.time()
m_linreg = LinearRegression().fit(X_train[linreg_train_cols].loc[train_inds], y_train.loc[train_inds])
t_linreg = time.time()-t0

#Fit Random Forest
t0 = time.time()
m_rf = RandomForestRegressor(n_estimators=100,max_depth=5, random_state=26, n_jobs=-1).fit(X_train.loc[train_inds], y_train.loc[train_inds])
t_rf = time.time()-t0/<code>

模型預測

值得注意的是，在訓練階段，我們使用了7 天前的 demand 數據以及 7 天前至 14 天之前的 demand 移動平均值數據作為特徵數據。但是在預測階段，是沒有 demand 數據的。因此這裡需要藉助滑動窗口，sliding window，的概念，也就是每次計算一個預測數據。 為了計算移動平均值數據，設置滑動窗口長度為 15。

通過滑動窗口方法預測未知數據的具體代碼如下：

<code>fday = datetime(2016,3, 28) 
max_lags = 15
for tdelta in range(0, 28):
    day = fday + timedelta(days=tdelta)
    tst = test[(test.date >= day - timedelta(days=max_lags)) & (test.date <= day)].copy()
    tst = feat_eng(tst)
    tst_lgb = tst.loc[tst.date == day , lgb_train_cols].copy()
    test.loc[test.date == day, "preds_LightGB"] = m_lgb.predict(tst_lgb)
    tst_rf = tst.loc[tst.date == day , lgb_train_cols].copy()
    tst_rf = tst_rf.fillna(0) 
    test.loc[test.date == day, "preds_RandomForest"] = m_rf.predict(tst_rf)
    
    tst_linreg = tst.loc[tst.date == day , linreg_train_cols].copy()
    tst_linreg = tst_linreg.fillna(0)    
    test.loc[test.date == day, "preds_LinearReg"] = m_linreg.predict(tst_linreg) 
    
test_final = test.loc[test.date >= fday]/<code>

Light Gradient Boosting

<code>model_name='LightGB'
predictions[model_name] = test_final["preds_"+model_name]

#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test_final.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test_final['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t_lgb, 'RMSE':score},ignore_index=True)/<code>

Light Gradient Boosting 預測結果

Random Forest

<code>model_name='RandomForest'
predictions[model_name] = test_final["preds_"+model_name]

#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test_final.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test_final['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t_lgb, 'RMSE':score},ignore_index=True)/<code>

Random Forest 預測結果

Linear Regression

<code>model_name='LinearReg'
predictions[model_name] = test_final["preds_"+model_name]

#visualize
fig, ax = plt.subplots(figsize=(25,4))
train[-28:].plot(x='date',y='demand',label='Train',ax=ax)
test_final.plot(x='date',y='demand',label='Test',ax=ax);
predictions.plot(x='date',y=model_name,label=model_name,ax=ax);
#evaluate
score = np.sqrt(mean_squared_error(predictions[model_name].values, test_final['demand']))
print('RMSE for {}: {:.4f}'.format(model_name,score))

stats = stats.append({'Model Name':model_name, 'Execution Time':t_linreg, 'RMSE':score},ignore_index=True)/<code>

Linear Regression 預測結果

以上就是所有的預測方法及過程。各個方法的運算時間及結果誤差如下：

<code>stats.sort_values(by='RMSE')
stats.plot(kind='bar',x='Model Name', y='RMSE', figsize=(12,6), title="Model RMSE Comparison - Lower is better");/<code>

各個方法的運算時間及結果誤差對比

各個方法的結果誤差對比

可以看出，傳統預測方法的性能相較於機器學習預測方法較差。

但是這個結論並不是絕對的。方法的準確度取決於不同的問題背景。 機器學習方法依賴於特徵數據。如果我們只有時間序列數據，那麼特徵數據較為缺乏，我們可以基於原始數據創建特徵數據，如滯後數據、移動平均數據等。因此機器學習方法要呈現更好地預測結果，特徵工程至關重要。在機器學習領域，某種程度上，數據才是起決定作用，而不是模型或者算法。

閱讀更多 deephub 的文章

關鍵字: 預測 社交網絡 機器學習