[譯] 時間序列分析、可視化、和使用 LSTM 預測

• 原文地址：Time Series Analysis, Visualization & Forecasting with LSTM
• 原文做者：Susan Li
• 譯文出自：掘金翻譯計劃
• 本文永久連接：github.com/xitu/gold-m…
• 譯者：Minghao23
• 校對者：Xuyuey，TrWestdoor

統計正態性檢驗，平穩性 Dickey-Fuller 檢驗，長短時間記憶網絡

標題已經闡述了一切。

閒話少說，讓咱們直接開始吧！

數據

該數據是在近四年的時間裏對一個家庭以一分鐘採樣率測量的電力消耗，能夠在這裏下載。

數據包括不一樣的電量值和一些分表的數值。然而，咱們只關注 Global_active_power 這個變量。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
pd.set_option('display.float_format', lambda x: '%.4f' % x)
import seaborn as sns
sns.set_context("paper", font_scale=1.3)
sns.set_style('white')
import warnings
warnings.filterwarnings('ignore')
from time import time
import matplotlib.ticker as tkr
from scipy import stats
from statsmodels.tsa.stattools import adfuller
from sklearn import preprocessing
from statsmodels.tsa.stattools import pacf
%matplotlib inline

import math
import keras
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.layers import Dropout
from keras.layers import *
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from keras.callbacks import EarlyStopping

df=pd.read_csv('household_power_consumption.txt', delimiter=';')
print('Number of rows and columns:', df.shape)
df.head(5)
複製代碼

如下數據預處理和特徵工程步驟須要完成：

• 將日期和時間合併到同一列，並轉換爲 datetime 類型。
• 將 Global_active_power 轉換爲數值型，並移除缺失值（1.2%）。
• 建立年、季度、月和日的特徵。
• 建立周的特徵，「0」表示週末，「1」表示工做日。

df['date_time'] = pd.to_datetime(df['Date'] + ' ' + df['Time'])
df['Global_active_power'] = pd.to_numeric(df['Global_active_power'], errors='coerce')
df = df.dropna(subset=['Global_active_power'])
df['date_time']=pd.to_datetime(df['date_time'])
df['year'] = df['date_time'].apply(lambda x: x.year)
df['quarter'] = df['date_time'].apply(lambda x: x.quarter)
df['month'] = df['date_time'].apply(lambda x: x.month)
df['day'] = df['date_time'].apply(lambda x: x.day)
df=df.loc[:,['date_time','Global_active_power', 'year','quarter','month','day']]
df.sort_values('date_time', inplace=True, ascending=True)
df = df.reset_index(drop=True)
df["weekday"]=df.apply(lambda row: row["date_time"].weekday(),axis=1)
df["weekday"] = (df["weekday"] < 5).astype(int)

print('Number of rows and columns after removing missing values:', df.shape)
print('The time series starts from: ', df.date_time.min())
print('The time series ends on: ', df.date_time.max())
複製代碼

移除缺失值以後，數據包括從 2006 年 12 月到 2010 年 11 月（47 個月）共 2,049,280 個測量值。

初始數據包括多個變量。這裏咱們只會關注一個單獨的變量：房屋的 Global_active_power 歷史記錄，也就是整個房屋平均每分鐘消耗的有功功率，單位是千瓦。

統計正態性檢驗

有一些統計測試方法能夠用來量化咱們的數據是否看起來像高斯分佈採樣。咱們將會使用 D’Agostino’s K² 檢驗。

在 SciPy 對這個檢驗的實現中，咱們對 p 值作出以下解釋。

• p <= alpha：拒絕 H0，非正態。
• p > alpha：不拒絕 H0，正態。

stat, p = stats.normaltest(df.Global_active_power)
print('Statistics=%.3f, p=%.3f' % (stat, p))
alpha = 0.05
if p > alpha:
    print('Data looks Gaussian (fail to reject H0)')
else:
    print('Data does not look Gaussian (reject H0)')
複製代碼

同時咱們也會計算峯度和偏度，以肯定數據分佈是否偏離正態分佈。

sns.distplot(df.Global_active_power);
print( 'Kurtosis of normal distribution: {}'.format(stats.kurtosis(df.Global_active_power)))
print( 'Skewness of normal distribution: {}'.format(stats.skew(df.Global_active_power)))
複製代碼

峯度：描述分佈的尾重

正態分佈的峯度接近於 0。若是峯度大於 0，則分佈尾部較重。若是峯度小於 0，則分佈尾部較輕。咱們計算出的峯度是大於 0 的。

偏度: 度量分佈的不對稱性

若是偏度介於 -0.5 和 0.5 之間，則數據是基本對稱的。若是偏度介於 -1 和 -0.5 之間或者 0.5 和 1 之間，則數據是稍微偏斜的。若是偏度小於 -1 或大於 1， 則數據是高度偏斜的。咱們計算出的偏度是大於 1 的。

第一個時間序列圖像

df1=df.loc[:,['date_time','Global_active_power']]
df1.set_index('date_time',inplace=True)
df1.plot(figsize=(12,5))
plt.ylabel('Global active power')
plt.legend().set_visible(False)
plt.tight_layout()
plt.title('Global Active Power Time Series')
sns.despine(top=True)
plt.show();
複製代碼

很明顯，這個圖像並非咱們想要的。不要這麼作。

年度和季度整體有功功率箱形圖對比

plt.figure(figsize=(14,5))
plt.subplot(1,2,1)
plt.subplots_adjust(wspace=0.2)
sns.boxplot(x="year", y="Global_active_power", data=df)
plt.xlabel('year')
plt.title('Box plot of Yearly Global Active Power')
sns.despine(left=True)
plt.tight_layout()

plt.subplot(1,2,2)
sns.boxplot(x="quarter", y="Global_active_power", data=df)
plt.xlabel('quarter')
plt.title('Box plot of Quarterly Global Active Power')
sns.despine(left=True)
plt.tight_layout();
複製代碼

當並排比較每一年的箱形圖時，咱們注意到 2006 年的整體有功功率的中位數相比於其餘年份高不少。其實這裏會有一點誤導。若是你還記得，咱們只有 2006 年 12 月的數據。而很明顯 12 月是一個家庭電力消耗的高峯月。

季度整體有功功率的中位數就比較符合預期，第1、四季度（冬季）較高，第三季度（夏季）最低。

整體有功功率分佈

plt.figure(figsize=(14,6))
plt.subplot(1,2,1)
df['Global_active_power'].hist(bins=50)
plt.title('Global Active Power Distribution')

plt.subplot(1,2,2)
stats.probplot(df['Global_active_power'], plot=plt);
df1.describe().T
複製代碼

正態機率圖也顯示這個數據與正態分佈偏離很大。

按照天、周、月、季度和年從新抽樣平均整體有功功率

fig = plt.figure(figsize=(18,16))
fig.subplots_adjust(hspace=.4)
ax1 = fig.add_subplot(5,1,1)
ax1.plot(df1['Global_active_power'].resample('D').mean(),linewidth=1)
ax1.set_title('Mean Global active power resampled over day')
ax1.tick_params(axis='both', which='major')

ax2 = fig.add_subplot(5,1,2, sharex=ax1)
ax2.plot(df1['Global_active_power'].resample('W').mean(),linewidth=1)
ax2.set_title('Mean Global active power resampled over week')
ax2.tick_params(axis='both', which='major')

ax3 = fig.add_subplot(5,1,3, sharex=ax1)
ax3.plot(df1['Global_active_power'].resample('M').mean(),linewidth=1)
ax3.set_title('Mean Global active power resampled over month')
ax3.tick_params(axis='both', which='major')

ax4  = fig.add_subplot(5,1,4, sharex=ax1)
ax4.plot(df1['Global_active_power'].resample('Q').mean(),linewidth=1)
ax4.set_title('Mean Global active power resampled over quarter')
ax4.tick_params(axis='both', which='major')

ax5  = fig.add_subplot(5,1,5, sharex=ax1)
ax5.plot(df1['Global_active_power'].resample('A').mean(),linewidth=1)
ax5.set_title('Mean Global active power resampled over year')
ax5.tick_params(axis='both', which='major');
複製代碼

一般來講，咱們的時間序列不會存在上升或降低的趨勢。最高的平均耗電量彷佛是在 2007 年以前，實際上這是由於咱們在 2007 年只有 12 月的數據（譯者注：原文有誤，應該是隻有 2006 年 12 月的數據），而那個月是用電高峯月。也就是說，若是咱們逐年比較，這個序列其實較爲平穩。

繪製整體有功功率均值圖，並以年、季、月和天分組

plt.figure(figsize=(14,8))
plt.subplot(2,2,1)
df.groupby('year').Global_active_power.agg('mean').plot()
plt.xlabel('')
plt.title('Mean Global active power by Year')

plt.subplot(2,2,2)
df.groupby('quarter').Global_active_power.agg('mean').plot()
plt.xlabel('')
plt.title('Mean Global active power by Quarter')

plt.subplot(2,2,3)
df.groupby('month').Global_active_power.agg('mean').plot()
plt.xlabel('')
plt.title('Mean Global active power by Month')

plt.subplot(2,2,4)
df.groupby('day').Global_active_power.agg('mean').plot()
plt.xlabel('')
plt.title('Mean Global active power by Day');
複製代碼

以上的圖像證明了咱們以前的發現。以年爲單位，序列較爲平穩。以季度爲單位，最低的平均耗電量處於第三季度。以月爲單位，最低的平均耗電量處於七月和八月。以天爲單位，最低的平均耗電量大約在每個月的 8 號（不知道爲何）。

每一年的整體有功功率

這一次，咱們移除 2006 年。

pd.pivot_table(df.loc[df['year'] != 2006], values = "Global_active_power",
               columns = "year", index = "month").plot(subplots = True, figsize=(12, 12), layout=(3, 5), sharey=True);
複製代碼

從 2007 年到 2010 年，每一年的模式都很類似。

工做日和週末的整體有功功率對比

dic={0:'Weekend',1:'Weekday'}
df['Day'] = df.weekday.map(dic)

a=plt.figure(figsize=(9,4))
plt1=sns.boxplot('year','Global_active_power',hue='Day',width=0.6,fliersize=3,
                    data=df)
a.legend(loc='upper center', bbox_to_anchor=(0.5, 1.00), shadow=True, ncol=2)
sns.despine(left=True, bottom=True)
plt.xlabel('')
plt.tight_layout()
plt.legend().set_visible(False);
複製代碼

在 2010 年之前，工做日的整體有功功率的中位數要比周末低一些。在 2010 年，它們徹底相等。

工做日和週末的整體有功功率對比的因素圖

plt1=sns.factorplot('year','Global_active_power',hue='Day',
                    data=df, size=4, aspect=1.5, legend=False)
plt.title('Factor Plot of Global active power by Weekend/Weekday')
plt.tight_layout()
sns.despine(left=True, bottom=True)
plt.legend(loc='upper right');
複製代碼

以年爲單位，工做日和週末都遵循一樣的模式。

原則上，當使用 LSTM 時，咱們不須要去檢驗或修正平穩性。然而，若是數據是平穩的，它會幫助模型提升性能，使神經網絡更容易學習。

平穩性

在統計學中，Dickey–Fuller test 檢驗了一個零假設，即單位根存在於自迴歸模型中。備擇假設依據使用的檢驗方法的不一樣而不一樣，可是一般爲平穩性或趨勢平穩性。

平穩序列的均值和方差一直是常數。時間序列在滑動窗口下的均值和標準差不隨時間變化。

Dickey-Fuller 檢驗

零檢驗（H0）：代表時間序列有一個單位根，意味着它是非平穩的。它包含一些和時間相關的成分。

備擇檢驗（H1）：代表時間序列不存在單位根，意味着它是平穩的。它不包含和時間相關的成分。

p-value > 0.05：接受零檢驗（H0），數據有單位根且是非平穩的。

p-value <= 0.05：拒絕零檢驗（H0），數據沒有單位根且是平穩的。

df2=df1.resample('D', how=np.mean)

def test_stationarity(timeseries):
    rolmean = timeseries.rolling(window=30).mean()
    rolstd = timeseries.rolling(window=30).std()

    plt.figure(figsize=(14,5))
    sns.despine(left=True)
    orig = plt.plot(timeseries, color='blue',label='Original')
    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
    std = plt.plot(rolstd, color='black', label = 'Rolling Std')

    plt.legend(loc='best'); plt.title('Rolling Mean & Standard Deviation')
    plt.show()

    print ('<Results of Dickey-Fuller Test>')
    dftest = adfuller(timeseries, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4],
                         index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
        dfoutput['Critical Value (%s)'%key] = value
    print(dfoutput)
test_stationarity(df2.Global_active_power.dropna())
複製代碼

從以上結論可得，咱們會拒絕零檢驗 H0，由於數據沒有單位根且是平穩的。

LSTM

咱們的任務是根據一個家庭兩百萬分鐘的耗電量歷史記錄，對這個時間序列作預測。咱們將使用一個多層的 LSTM 遞歸神經網絡來預測時間序列的最後一個值。

若是你想縮減計算時間，並快速得到結果來檢驗模型，你能夠對數據以小時爲單位從新採樣。在本文的實驗中我會維持原單位爲分鐘。

在構建 LSTM 模型以前，須要進行下列數據預處理和特徵工程的工做。

• 建立數據集，保證全部的數據的類型都是 float。
• 特徵標準化。
• 分割訓練集和測試集。
• 將數值數組轉換爲數據集矩陣。
• 將維度轉化爲 X=t 和 Y=t+1。
• 將輸入維度轉化爲三維 (num_samples, num_timesteps, num_features)。

dataset = df.Global_active_power.values #numpy.ndarray
dataset = dataset.astype('float32')
dataset = np.reshape(dataset, (-1, 1))
scaler = MinMaxScaler(feature_range=(0, 1))
dataset = scaler.fit_transform(dataset)
train_size = int(len(dataset) * 0.80)
test_size = len(dataset) - train_size
train, test = dataset[0:train_size,:], dataset[train_size:len(dataset),:]

def create_dataset(dataset, look_back=1):
    X, Y = [], []
    for i in range(len(dataset)-look_back-1):
        a = dataset[i:(i+look_back), 0]
        X.append(a)
        Y.append(dataset[i + look_back, 0])
    return np.array(X), np.array(Y)

look_back = 30
X_train, Y_train = create_dataset(train, look_back)
X_test, Y_test = create_dataset(test, look_back)

# 將輸入維度轉化爲 [samples, time steps, features]
X_train = np.reshape(X_train, (X_train.shape[0], 1, X_train.shape[1]))
X_test = np.reshape(X_test, (X_test.shape[0], 1, X_test.shape[1]))
複製代碼

模型結構

• 定義 LSTM 模型，第一個隱藏層含有 100 個神經元，輸出層含有 1 個神經元，用於預測 Global_active_power。輸入的維度是一個包含 30 個特徵的時間步長。
• Dropout 20%。
• 使用均方差損失函數，和改進於隨機梯度降低的效率更高的 Adam。
• 模型將會進行 20 個 epochs 的訓練，每一個 batch 的大小爲 70。

model = Sequential()
model.add(LSTM(100, input_shape=(X_train.shape[1], X_train.shape[2])))
model.add(Dropout(0.2))
model.add(Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')

history = model.fit(X_train, Y_train, epochs=20, batch_size=70, validation_data=(X_test, Y_test),
                    callbacks=[EarlyStopping(monitor='val_loss', patience=10)], verbose=1, shuffle=False)

model.summary()
複製代碼

作出預測

train_predict = model.predict(X_train)
test_predict = model.predict(X_test)
# 預測值求逆
train_predict = scaler.inverse_transform(train_predict)
Y_train = scaler.inverse_transform([Y_train])
test_predict = scaler.inverse_transform(test_predict)
Y_test = scaler.inverse_transform([Y_test])

print('Train Mean Absolute Error:', mean_absolute_error(Y_train[0], train_predict[:,0]))
print('Train Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_train[0], train_predict[:,0])))
print('Test Mean Absolute Error:', mean_absolute_error(Y_test[0], test_predict[:,0]))
print('Test Root Mean Squared Error:',np.sqrt(mean_squared_error(Y_test[0], test_predict[:,0])))
複製代碼

繪製模型損失

plt.figure(figsize=(8,4))
plt.plot(history.history['loss'], label='Train Loss')
plt.plot(history.history['val_loss'], label='Test Loss')
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epochs')
plt.legend(loc='upper right')
plt.show();
複製代碼

比較真實值和預測值

在個人結果中，每一個時間步長是 1 分鐘。若是你以前以小時從新採樣了數據，那麼在你的結果裏每一個時間步長是 1 小時。

我將會比較最近 200 分鐘的真實值和預測值。

aa=[x for x in range(200)]
plt.figure(figsize=(8,4))
plt.plot(aa, Y_test[0][:200], marker='.', label="actual")
plt.plot(aa, test_predict[:,0][:200], 'r', label="prediction")
# plt.tick_params(left=False, labelleft=True) # 移除 ticks
plt.tight_layout()
sns.despine(top=True)
plt.subplots_adjust(left=0.07)
plt.ylabel('Global_active_power', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show();
複製代碼

LSTMs 太神奇了！

Jupyter notebook 能夠在 Github 中找到。享受這一週餘下的時光吧！

參考： Multivariate Time Series Forecasting with LSTMs in Keras

若是發現譯文存在錯誤或其餘須要改進的地方，歡迎到 掘金翻譯計劃 對譯文進行修改並 PR，也可得到相應獎勵積分。文章開頭的 本文永久連接 即爲本文在 GitHub 上的 MarkDown 連接。

---

掘金翻譯計劃 是一個翻譯優質互聯網技術文章的社區，文章來源爲 掘金 上的英文分享文章。內容覆蓋 Android、iOS、前端、後端、區塊鏈、產品、設計、人工智能等領域，想要查看更多優質譯文請持續關注 掘金翻譯計劃、官方微博、知乎專欄。