本文以Kaggle上一個能源需求預測的案例為基礎，實戰示範指數平滑（Holt-Winters）方法應用全流程。調用模型雖簡單，但是入模前的資料分析、處理，模型參數的優化、效果的分析亦尤為重要，能夠分析全面産出最優的方案又顯得不那麼簡單。本文期望通過一個簡單場景，回顧ML應用的基本流程。
導入子產品
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.holtwinters import ExponentialSmoothing
讀取資料
來自PJM超過10年的每小時能源消耗資料，機關：兆瓦。
df = pd.read_csv('data/PJME_hourly.csv')
df['Datetime'] = pd.to_datetime(df['Datetime'])
df.sort_values(by=['Datetime'], axis=0, ascending=True, inplace=True)
df.reset_index(inplace=True, drop=True)
print(df)
可視化
df.set_index('Datetime').plot(style='r', figsize=(20,4), legend=False)
plt.show()
檢查資料
從上圖資料基本已經能看出資料時間範圍、資料粒度、資料大體特征（周期型）、直覺上是否有異常點（2012～2014貌似有一個）。
>> 那就先來看看資料是否有缺失值 -> 無。
print(df[pd.isnull(df['PJME_MW'])])
>> 再來看看是否有異常值，結合資料可視化時的序列曲線圖感覺2012～2014可能有異常值。進一步探索發現異常偏低資料集中在12年的10月30日當天。
print(df[df['PJME_MW'] <= 15000])
>> 但同往年趨勢特征對比整體水準偏低但并無明顯異常，是以決定對此不再做異常處理。
ts1 = df[(df['Datetime'] >= '2012-10-29 01:00:00') & (df['Datetime'] <= '2012-10-31 23:00:00')]['PJME_MW']
ts2 = df[(df['Datetime'] >= '2011-10-29 01:00:00') & (df['Datetime'] <= '2011-10-31 23:00:00')]['PJME_MW']
plt.figure(figsize=(20,4))
ts1.reset_index(drop=True).plot(label='2012', legend=True)
ts2.reset_index(drop=True).plot(label='2011', legend=True)
plt.show()
>> 下面接着再來看看時間粒度（讀取資料時能看出來應該是小時級粒度）、是否等頻（有無日期缺失、有無日期重複）。發現資料日期有重複且有日期缺失情況。（下圖中黃色字型均值為1表明1小時1條）。
df['SPAN'] = df['Datetime'].diff().astype('timedelta64[h]').astype(float)
print(df.head())
print(df['SPAN'].median(), df['SPAN'].min(), df['SPAN'].max())
>> 既然資料日期有重複、有遺漏，不免好奇哪些重複、哪些遺漏，或者說重複&遺漏下資料的表現是怎樣的。
# 日期重複資料
print(pd.merge(df, df[df['SPAN'] == 0]['Datetime'], on='Datetime'))
# 漏報前後資料
index = df[df['SPAN'] == 2].index[0]
print(df.loc[index-3: index+3])
# 漏報次數
print(df[df['SPAN'] == 2].shape[0]
清洗資料
根據以上檢查，發現資料相對還是比較幹淨的，但有2點問題：
1.日期有重複 -> 以上檢查中可以看到資料并不是完全重複，猜測可能偶有一分鐘内多次上傳資料，此時為了資料粒度統一，舍棄較早的一條。
2.資料有漏報 -> 因傳感器/網絡等設施問題難免有資料漏報的情況，目前資料漏報間隔較短且漏報次數較少，可以使用插值的方法補充漏報資料。
# 去重
df.drop_duplicates(subset='Datetime', keep='last', inplace=True)
# 插值
df.set_index('Datetime', inplace=True)
date_range = pd.date_range(start=df.index.min(), end=df.index.max(), freq='H')
df = df.reindex(date_range)
df['PJME_MW'].interpolate(method='linear', inplace=True)
print(df)
# 删除臨時列
df.drop('SPAN', axis=1, inplace=True)
拆分訓練集/驗證集
需求為預測未來1年的資料，首先一個問題是不管用什麼模型都至少要有一份資料能夠驗證模型效果。是以首先要拆分訓練集和驗證集，如果要更精确的驗證模型效果還需拆分測試集。時間序列中一般根據時間拆分，較早的一段時間資料用于訓練模型，最近的一段時間用于驗證/測試模型效果。
比如我們要用過去5年的資料訓練模型，可以留出最近1年的資料作為驗證集，再往前推5年的資料做訓練集。
train_years = 5
val_years = 1
last_date = df.index[-1]
val_start_date = last_date - pd.DateOffset(years=val_years)
train_start_date = val_start_date - pd.DateOffset(months=train_years)
print(train_start_date, val_start_date, last_date)
df_train = df.loc[(df.index >= train_start_date) & (df.index < val_start_date)]
df_val = df.loc[df.index >= val_start_date]
訓練模型
模型選擇
從資料的可視化圖中能直覺的看出來，目前資料為無趨勢的周期型時間序列資料，周期大小為1年。本次限定使用指數平滑方法，有周期故使用三次指數平滑（Holt-Winters）方法，具體使用季節性加法模型還是乘法模型可進行實驗根據評估名額選擇。
模型訓練
>> 評估名額：MAPE
def mape(y_true, y_pred):
y_true, y_pred = np.array(y_true), np.array(y_pred)
pe = (y_true - y_pred) / y_true
ape = np.abs(pe)
mape = np.mean(ape)
return f'{mape*100:.2f}%'
>> 訓練模型
周期大小為24*365（一年，每天24個資料點），無趨勢不指定trend，有季節性暫指定為加法模型。
模型訓練時考慮資料量較大，設定模型參數初始值不進行暴力搜尋，優化算法使用‘L-BFGS-B’，目的為提高計算速度，不然實在太慢。
periods = 24*365
fit = ExponentialSmoothing(df_train['PJME_MW'],
seasonal_periods=periods,
seasonal='add').fit(method='L-BFGS-B', use_brute=False)
print(fit)
fcast = fit.forecast(len(df_val))
mape_hw = mape(y_true=df_val['PJME_MW'], y_pred=fcast)
print(f'Our Holt-Winter model has a mean average percentage error of {mape_hw}')
>> 參數調優
HW中的參數使用‘L-BFGS-B’優化算法自動優化，對比季節性加法模型，乘法模型效果較好MAPE為11.75%。
fit = ExponentialSmoothing(df_train['PJME_MW'],
seasonal_periods=periods,
seasonal='mul').fit(method='L-BFGS-B', use_brute=False)
fcast = fit.forecast(len(df_val))
mape_hw = mape(y_true=df_val['PJME_MW'], y_pred=fcast)
print(f'Our Holt-Winter model has a mean average percentage error of {mape_hw1}')
模型分析
通過進一步分析，認為使用7年資料訓練模型效果較好，此時在近1年上的驗證集上MAPE為9.7%。
分析内容有以下兩點：
(1). 分析訓練資料的大小對模型效果的影響
(2). 分析模型在不同年份上的表現看模型是否穩定
def test_train_size(df):
val_years = 1
periods = 24*365
for train_years in range(1,10):
# 拆分日期
last_date = df.index[-1]
val_start_date = last_date - pd.DateOffset(years=val_years)
train_start_date = val_start_date - pd.DateOffset(years=train_years)
# 訓練集/驗證集
df_train = df.loc[(df.index >= train_start_date) & (df.index < val_start_date)]
df_val = df.loc[df.index >= val_start_date]
# 模型訓練
fit = ExponentialSmoothing(df_train['PJME_MW'],
seasonal_periods=periods,
seasonal='mul').fit(method='L-BFGS-B', use_brute=False)
# 模型預測
fcast = fit.forecast(len(df_val))
# 模型評估
mape_hw = mape(y_true=df_val['PJME_MW'], y_pred=fcast)
# 列印結果
if train_years == 1:
val_date_desc = 'val_date_range : {}-{}'.format(val_start_date.strftime('%Y%m%d'), last_date.strftime('%Y%m%d'))
print(val_date_desc)
train_date_desc = 'train_date_range: {}-{}'.format(train_start_date.strftime('%Y%m%d'), val_start_date.strftime('%Y%m%d'))
print(train_date_desc, 'train_years:', train_years, 'mape:', mape_hw)
# 分析不同年數訓練資料在18年（20170803-20180803）驗證集上的表現
test_train_size(df)
# 分析不同年數訓練資料在17年（20160803-20170803）驗證集上的表現
last_date = df.index[-1] - pd.DateOffset(years=1)
test_train_size(df.loc[df.index<=last_date])
模型預測
# 起始日期
last_date = df.index[-1]
train_start_date = last_date - pd.DateOffset(years=7)
# 訓練集
df_train = df.loc[df.index >= train_start_date]
# 模型訓練
fit = ExponentialSmoothing(df_train['PJME_MW'],
seasonal_periods=periods,
seasonal='mul').fit(method='L-BFGS-B', use_brute=False)
# 模型預測
fcast = fit.forecast(24*365)
date_range = pd.date_range(start=df_train.index.min(), end=fcast.index.max(), freq='H')
df_train = df_train.reindex(date_range)
df_train['predict'] = fcast
print(fcast.index.min(), fcast.index.max())
# 可視化
plt.figure(figsize=(20,4))
df_train['PJME_MW'].plot()
df_train['predict'].plot()
plt.xlabel('Date & Time')
plt.ylabel('Energy Demand [MW]')
plt.title('Holt-Winter Forecast of Hourly Energy Demand')
plt.show()