[Перевод] Анализ временных рядов, визуализация и прогнозирование с помощью LSTM

машинное обучение искусственный интеллект Программа перевода самородков Python задняя часть
[Перевод] Анализ временных рядов, визуализация и прогнозирование с помощью LSTM

Статистический критерий нормальности, критерий стационарности Дики-Фуллера, сети с долговременной кратковременной памятью

Название говорит само за себя.

Без лишних слов, давайте начнем!

данные

Данные представляют собой потребление электроэнергии, измеренное при одноминутной частоте выборки для домохозяйства в течение почти четырех лет, и их можно найти по адресуздесьскачать.

Данные включают в себя различные значения мощности и некоторые субметровые значения. Однако мы сосредоточимся только на переменной 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)

Table 1

Необходимо выполнить следующие этапы предварительной обработки данных и разработки признаков:

  • Объедините дату и время в один и тот же столбец и преобразуйте их в тип 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())

После удаления пропущенных значений данные включали 2 049 280 измерений с декабря 2006 г. по ноябрь 2010 г. (47 месяцев).

Исходные данные включают несколько переменных. Здесь мы сосредоточимся только на одной переменной: истории global_active_power дома, которая представляет собой среднюю активную мощность, потребляемую всем домом в минуту в киловаттах.

Статистический тест на нормальность

Есть некоторые статистические тесты, которые мы можем использовать, чтобы количественно определить, выглядят ли наши данные как выборка по Гауссу. мы будем использоватьK²-тест Д’Агостино.

существуетSciPyПри реализации этого теста мы интерпретируем значение p следующим образом.

  • p
  • р > альфа: 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)))

图 1

эксцесс: описывает хвостовой вес распределения

Эксцесс нормального распределения близок к 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();

图 2

Очевидно, что это изображение не то, что нам нужно. Не делай этого.

Сравнение годовой и квартальной общей диаграммы Active Power Box

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

图 3

Сравнивая диаграммы для каждого года рядом, мы замечаем, что медиана общей активной мощности в 2006 г. намного выше, чем в другие годы. На самом деле здесь есть некоторое заблуждение. Если вы помните, у нас есть данные только за декабрь 2006 года. И очевидно, что декабрь — это пиковый месяц потребления электроэнергии в домашних хозяйствах.

Медианная квартальная общая активная мощность больше соответствует ожиданиям, выше в первом и четвертом кварталах (зима) и самая низкая в третьем квартале (лето).

Общее распределение активной мощности

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

图 4

График нормальной вероятности также показывает, что эти данные значительно отклоняются от нормального распределения.

Повторная выборка средней общей активной мощности по дням, неделям, месяцам, кварталам и годам

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

图 5

Вообще говоря, наш временной ряд не будет иметь восходящей или нисходящей тенденции. Самое высокое среднее потребление электроэнергии, по-видимому, было до 2007 года, на самом деле это потому, что у нас есть данные только за декабрь 2007 года, а этот месяц является пиковым месяцем потребления электроэнергии. То есть, если сравнивать по годам, то серия на самом деле относительно стабильна.

Постройте среднюю общую активную мощность, сгруппированную по годам, кварталам, месяцам и дням.

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

图 6

Изображения выше подтверждают наши предыдущие выводы. По годам ряд относительно стабилен. В квартальном исчислении самое низкое среднее энергопотребление приходится на третий квартал. В месячном исчислении самое низкое среднее потребление электроэнергии приходится на июль и август. В днях самое низкое среднее энергопотребление приходится на 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);

图 7

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

图 8

До 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');

图 9

Ежегодно и будни, и выходные следуют одной и той же схеме.

В принципе, при использованииLSTM, нам не нужно проверять или исправлятьСтационарность. Однако, если данные стационарны, это поможет модели повысить производительность и упростить обучение нейронной сети.

Стационарность

В статистике,Тест Дики-ФуллераПроверялась нулевая гипотеза о существовании единичного корня в авторегрессионной модели. Альтернативные гипотезы варьируются в зависимости от используемого метода тестирования, но обычноСтационарностьилистационарность тренда.

Среднее значение и дисперсия стационарного ряда всегда постоянны. Среднее значение и стандартное отклонение временного ряда в скользящем окне не меняются со временем.

Тест Дики-Фуллера

нулевой тест(H0): указывает на то, что временной ряд имеет единичный корень, то есть он нестационарен. Он содержит некоторые компоненты, связанные со временем.

Альтернативный тест(H1): указывает на то, что временной ряд не имеет единичного корня, то есть является стационарным. Он не содержит компонентов, связанных со временем.

p-значение> 0,05: критерий нуля (H0) принят, данные имеют единичный корень и нестационарны.

p-значение 

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

图 10

Из приведенного выше вывода мы отклоним нулевой критерий H0, поскольку данные не имеют единичного корня и являются стационарными.

LSTM

Наша задача состоит в том, чтобы сделать прогнозы по этому временному ряду на основе 2 миллионов минут истории потребления электроэнергии домохозяйством. Мы будем использовать многослойную рекуррентную нейронную сеть 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 признаков.
  • Отсев 20%.
  • Используйте среднеквадратичную функцию потерь и более эффективный Адам, который улучшает стохастический градиентный спуск.
  • Модель будет обучаться в течение 20 эпох, каждая с размером пакета 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();

图 11

Сравните фактические и прогнозируемые значения

В моих результатах каждый временной шаг равен 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();

图 12

LSTM потрясающие!

Jupyter notebookдопустимыйGithubнайти в. Наслаждайтесь остальной частью недели!

Ссылаться на:Multivariate Time Series Forecasting with LSTMs in Keras

Если вы обнаружите ошибки в переводе или в других областях, требующих доработки, добро пожаловать наПрограмма перевода самородковВы также можете получить соответствующие бонусные баллы за доработку перевода и PR. начало статьиПостоянная ссылка на эту статьюЭто ссылка MarkDown этой статьи на GitHub.

Программа перевода самородковэто сообщество, которое переводит высококачественные технические статьи из Интернета сНаггетсДелитесь статьями на английском языке на . Охват контентаAndroid,iOS,внешний интерфейс,задняя часть,блокчейн,продукт,дизайн,искусственный интеллектЕсли вы хотите видеть более качественные переводы, пожалуйста, продолжайте обращать вниманиеПрограмма перевода самородков,официальный Вейбо,Знай колонку.

Категории