该数据是在近四年的时间里对一个家庭以一分钟采样率测量的电力消耗，可以在这里下载。
数据包括不同的电量值和一些分表的数值。然而，我们只关注 Global_active_power 这个变量。
import numpy as npimport matplotlib.pyplot as pltimport pandas as pdpd.set_option('display.float_format', lambda x: '%.4f' % x)import seaborn as snssns.set_context("paper", font_scale=1.3)sns.set_style('white')import warningswarnings.filterwarnings('ignore')from time import timeimport matplotlib.ticker as tkrfrom scipy import statsfrom statsmodels.tsa.stattools import adfullerfrom sklearn import preprocessingfrom statsmodels.tsa.stattools import pacf%matplotlib inlineimport mathimport kerasfrom keras.models import Sequentialfrom keras.layers import Densefrom keras.layers import LSTMfrom keras.layers import Dropoutfrom keras.layers import *from sklearn.preprocessing import MinMaxScalerfrom sklearn.metrics import mean_squared_errorfrom sklearn.metrics import mean_absolute_errorfrom keras.callbacks import EarlyStoppingdf=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.05if 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 月是一个家庭电力消耗的高峰月。
季度总体有功功率的中位数就比较符合预期，第一、四季度（冬季）较高，第三季度（夏季）最低。
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 检验了一个零假设，即单位根存在于自回归模型中。备择假设依据使用的检验方法的不同而不同，但是通常为平稳性或趋势平稳性。
平稳序列的均值和方差一直是常数。时间序列在滑动窗口下的均值和标准差不随时间变化。
零检验（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 模型之前，需要进行下列数据预处理和特征工程的工作。
创建数据集，保证所有的数据的类型都是 float。
特征标准化。
分割训练集和测试集。
将数值数组转换为数据集矩阵。
将维度转化为 X=t 和 Y=t+1。
将输入维度转化为三维 (num_samples, num_timesteps, num_features)。
dataset = df.Global_active_power.values #numpy.ndarraydataset = 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_sizetrain, 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 = 30X_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) # 移除 ticksplt.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 中找到。享受这一周余下的时光吧！
联系客服
