本文的出發點在于比較傳統預測方法和機器學習預測方法。

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

涉及到的方法有：

• 單指數平滑法
• 雙指數平滑法
• 三指數平滑法
• ARIMA
• SARIMA
• SARIMAX
• Light Gradient Boosting
• Random Forest
• Linear Regression

---

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

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'])

然后導入數據集：

data = pd.read_csv('data_for_tsa.csv') data['date'] = pd.to_datetime(data['date']) data.head()

數據集前五行數據

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

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

---

接下來進行數據集劃分

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

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);

時間序列數據

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

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

訓練及評價模型

單指數平滑方法

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

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

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

自相關圖

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

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

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)

單指數平滑方法預測結果

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

雙指數平滑方法

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

具體代碼如下：

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)

雙指數平滑方法預測結果

三指數平滑方法

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

具體代碼如下：

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)

三指數平滑方法預測結果

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

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 worksuppress_warnings=True, # we don't want convergence warningsstepwise=True) arima_model.summary()

auto_arima 的計算結果

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

#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)

ARIMA 預測結果

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

SARIMA

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

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

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 worksuppress_warnings=True, # we don't want convergence warningsstepwise=True) sarima_model.summary()

auto_arima 計算結果

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

#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)

SARIMA 預測結果

SARIMAX

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

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

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 worksuppress_warnings=True, # we don't want convergence warningsstepwise=True) sarimax_model.summary()

auto_arima 計算結果

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

#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)

SARIMAX 預測結果

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

機器學習

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

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

• 特征數據1：滯后數據。 選擇 7 天前的 demand 數據作為特征數據。
• 特征數據2：移動平均數據。選擇 7 天前至 14 天之前的 demand 移動平均值數據作為特征數據。
• 特征數據3：月銷售均值
• 特征數據4：每月銷售最大值
• 特征數據5：每月銷售最小值
• 特征數據6：每月銷售最大值與最小值的差值
• 特征數據7：每周銷售均值
• 特征數據8：每周銷售最大值
• 特征數據9：每周銷售中值

具體代碼如下：

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 dfdef per_timeframe_stats(df, col):#For each item compute its mean and other descriptive statistics for each month and dayofweek in the datasetmonths = 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 dfdef feat_eng(df):df = lags_windows(df)df = per_timeframe_stats(df,'demand')return df

準備數據：

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"]

模型擬合

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

#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

模型預測

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

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

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]

Light Gradient Boosting

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)

Light Gradient Boosting 預測結果

Random Forest

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)

Random Forest 預測結果

Linear Regression

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)

Linear Regression 預測結果

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

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

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

各個方法的結果誤差對比

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

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

最后，參加過這個比賽的應該都能看到我們的公眾號

總結

以上是生活随笔為你收集整理的KaggleM5 Forecasting：传统预测方法与机器学习预测方法对比的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。