Интуитивное использование сезонности для повышения точности модели
Добро пожаловать во вторую часть анализа временных рядов! В этом посте мы будем работать над моделированием данных временных рядов. Это продолжение моей предыдущей публикации о данных временных рядов.
В нашем предыдущем сообщении в блоге мы говорили о том, что такое данные временных рядов, как отформатировать такие данные, чтобы максимизировать их полезность, и как обрабатывать недостающие данные. Мы также узнали, как пересчитывать данные временных рядов по месяцам, неделям, годам и т. Д. И вычислять скользящие средние. Мы углубились в такие концепции, как тренды, сезонность, дифференциация первого порядка и автокорреляция. Если вы знакомы с большинством вещей, то все готово. Если вам нужно освежиться, вы можете выполнить быстрый поиск в Google по этим темам или прочитать мой предыдущий пост здесь.
Несколько слов перед тем, как мы начнем:
Существуют и другие, несомненно, лучшие пакеты, доступные для прогнозирования временных рядов, такие как ARIMA или собственное программное обеспечение Facebook Prophet. Однако эта статья была вдохновлена домашним заданием подруги, которое требовало, чтобы она использовала только Scikit, Numpy и Pandas (или столкнулась с мгновенной дисквалификацией!).
Давайте погрузимся в наш набор данных
Мы будем работать с общедоступным набором данных Open Power System Data. Вы можете скачать данные здесь. Он включает потребление электроэнергии, производство энергии ветра и производство солнечной энергии за 2006–2017 годы.
url='https://raw.githubusercontent.com/jenfly/opsd/master/opsd_germany_daily.csv' data = pd.read_csv(url,sep=",")
После установки столбца Date в качестве индекса наш набор данных выглядит так:
# to explicitly convert the date column to type DATETIME
data['Date'] = pd.to_datetime(data['Date'])
data = data.set_index('Date')
Определение задачи моделирования
Цели предсказания
Наша цель - предсказать Consumption (в идеале для будущих невидимых дат) на основе этого набора данных временного ряда.
Набор для обучения и тестирования
Мы будем использовать данные за 10 лет для обучения, то есть за 2006–2016 годы, и данные за прошлый год для тестирования, то есть за 2017 год.
Показатель производительности
Чтобы оценить, насколько хороша наша модель, мы будем использовать R-квадрат и среднеквадратичную ошибку (но будем печатать все соответствующие показатели, чтобы вы могли принять окончательное решение).
Вспомогательные функции
Чтобы распечатать все показатели производительности, относящиеся к задаче регрессии (например, MAE и R-квадрат), мы будем определять функцию regression_results.
import sklearn.metrics as metrics
def regression_results(y_true, y_pred):
# Regression metrics
explained_variance=metrics.explained_variance_score(y_true, y_pred)
mean_absolute_error=metrics.mean_absolute_error(y_true, y_pred)
mse=metrics.mean_squared_error(y_true, y_pred)
mean_squared_log_error=metrics.mean_squared_log_error(y_true, y_pred)
median_absolute_error=metrics.median_absolute_error(y_true, y_pred)
r2=metrics.r2_score(y_true, y_pred)
print('explained_variance: ', round(explained_variance,4))
print('mean_squared_log_error: ', round(mean_squared_log_error,4))
print('r2: ', round(r2,4))
print('MAE: ', round(mean_absolute_error,4))
print('MSE: ', round(mse,4))
print('RMSE: ', round(np.sqrt(mse),4))
Функциональная инженерия
В качестве основы мы выбираем упрощенную модель, которая прогнозирует сегодняшнюю стоимость потребления на основе
- вчерашнее значение потребления и;
- разница между вчерашним и позавчерашним значениями потребления.
# creating new dataframe from consumption column data_consumption = data[['Consumption']] # inserting new column with yesterday's consumption values data_consumption.loc[:,'Yesterday'] = data_consumption.loc[:,'Consumption'].shift() # inserting another column with difference between yesterday and day before yesterday's consumption values. data_consumption.loc[:,'Yesterday_Diff'] = data_consumption.loc[:,'Yesterday'].diff() # dropping NAs data_consumption = data_consumption.dropna()
Определение обучающих и тестовых наборов
X_train = data_consumption[:'2016'].drop(['Consumption'], axis = 1) y_train = data_consumption.loc[:'2016', 'Consumption'] X_test = data_consumption['2017'].drop(['Consumption'], axis = 1) y_test = data_consumption.loc['2017', 'Consumption']
Перекрестная проверка данных временных рядов
Во время собеседований по науке о данных часто возникает вопрос: Какой метод перекрестной проверки вы бы использовали для данных временных рядов?
У вас может возникнуть соблазн склониться к самой популярной кросс-валидации K-Fold (поверьте мне, до недавнего времени - не спрашивайте, сколько времени! - я не знал, что существуют другие методы CV, кроме K-fold). К сожалению, это был бы неправильный ответ. Причина в том, что он не принимает во внимание, что данные временных рядов имеют некоторый естественный порядок, а рандомизация при стандартной k-кратной перекрестной проверке не сохраняет этот порядок.
Лучшей альтернативой для перекрестной проверки данных временных рядов (чем K-кратное CV) является стратегия Forward Chaining.
В прямой цепочке, скажем, с 3 сгибами, наборы поездов и валидации выглядят так:
- кратность 1: обучение [1], проверка [2]
- кратность 2: обучение [1 2], проверка [3]
- кратность 3: обучение [1 2 3], проверка [4]
где 1, 2, 3, 4 представляют год. Таким образом, последовательные тренировочные наборы являются надмножествами предшествующих.
К счастью для нас, в sklearn есть возможность реализовать такое разделение теста поездов с использованием TimeSeriesSplit.
from sklearn.model_selection import TimeSeriesSplit
Функция TimeSerieSplit принимает в качестве входных данных количество разделений. Поскольку наши данные обучения имеют 11 уникальных лет (2006-2016), мы должны установить n_splits = 10. Таким образом, у нас есть аккуратные наборы для обучения и проверки:
- кратность 1: обучение [2006], валидация [2007]
- кратное 2: обучение [2006 2007], валидация [2008]
- тройка: обучение [2006 2007 2008], валидация [2009]
- 4 раза: обучение [2006 2007 2008 2009], валидация [2010]
- раз 5: обучение [2006 2007 2008 2009 2010], валидация [2011]
- раз 6: обучение [2006 2007 2008 2009 2010 2011], валидация [2012]
- раз 7: обучение [2006 2007 2008 2009 2010 2011 2012], проверка [2013]
- кратная 8: обучение [2006 2007 2008 2009 2010 2011 2012 2013], проверка [2014]
- раз 9: обучение [2006 2007 2008 2009 2010 2011 2012 2013 2014], проверка [2015]
- кратная 10: обучение [2006 2007 2008 2009 2010 2011 2012 2013 2014 2015], проверка [2016]
Алгоритмы выборочной проверки
# Spot Check Algorithms
models = []
models.append(('LR', LinearRegression()))
models.append(('NN', MLPRegressor(solver = 'lbfgs'))) #neural network
models.append(('KNN', KNeighborsRegressor()))
models.append(('RF', RandomForestRegressor(n_estimators = 10))) # Ensemble method - collection of many decision trees
models.append(('SVR', SVR(gamma='auto'))) # kernel = linear
# Evaluate each model in turn
results = []
names = []
for name, model in models:
# TimeSeries Cross validation
tscv = TimeSeriesSplit(n_splits=10)
cv_results = cross_val_score(model, X_train, y_train, cv=tscv, scoring='r2')
results.append(cv_results)
names.append(name)
print('%s: %f (%f)' % (name, cv_results.mean(), cv_results.std()))
# Compare Algorithms
plt.boxplot(results, labels=names)
plt.title('Algorithm Comparison')
plt.show()
И KNN, и RF работают одинаково хорошо. Но я лично предпочитаю RF, поскольку эта ансамблевая модель (объединяет несколько «индивидуальных» (разнообразных) моделей вместе и обеспечивает превосходную мощность прогнозирования.) Может работать практически сразу же, и это одна из причин, почему они очень популярны.
Гиперпараметры поиска в сетке
О необходимости поиска по сетке гиперпараметров я рассказывал в одной из своих предыдущих статей.
Оптимальная комбинация гиперпараметров максимизирует производительность модели, не приводя к проблеме высокой дисперсии (переобучения).
Код Python для выполнения поиска по сетке выглядит следующим образом:
from sklearn.model_selection import GridSearchCV
model = RandomForestRegressor()
param_search = {
'n_estimators': [20, 50, 100],
'max_features': ['auto', 'sqrt', 'log2'],
'max_depth' : [i for i in range(5,15)]
}
tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)
gsearch.fit(X_train, y_train)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_
Если вы заметили приведенный выше код, мы определили настраиваемый счетчик, установив scoring = rmse_score вместо использования одной из общих показателей оценки, определенных в sklearn. Мы определяем нашего индивидуального счетчика следующим образом:
from sklearn.metrics import make_scorer
def rmse(actual, predict):
predict = np.array(predict)
actual = np.array(actual)
distance = predict - actual
square_distance = distance ** 2
mean_square_distance = square_distance.mean()
score = np.sqrt(mean_square_distance)
return score
rmse_score = make_scorer(rmse, greater_is_better = False)
Проверка наилучшей производительности модели на тестовых данных
y_true = y_test.values y_pred = best_model.predict(X_test) regression_results(y_true, y_pred)
Для начала это неплохо. Посмотрим, сможем ли мы улучшить нашу модель.
Возврат при проектировании функций
До сих пор мы использовали значения на (t-1)th день для прогнозирования значений на t день. Теперь давайте также будем использовать значения от (t-2)days для прогнозирования потребления:
# creating copy of original dataframe data_consumption_2o = data_consumption.copy() # inserting column with yesterday-1 values data_consumption_2o['Yesterday-1'] = data_consumption_2o['Yesterday'].shift() # inserting column with difference in yesterday-1 and yesterday-2 values. data_consumption_2o['Yesterday-1_Diff'] = data_consumption_2o['Yesterday-1'].diff() # dropping NAs data_consumption_2o = data_consumption_2o.dropna()
Сброс поезда и тестового набора
X_train_2o = data_consumption_2o[:'2016'].drop(['Consumption'], axis = 1) y_train_2o = data_consumption_2o.loc[:'2016', 'Consumption'] X_test = data_consumption_2o['2017'].drop(['Consumption'], axis = 1) y_test = data_consumption_2o.loc['2017', 'Consumption']
Проверка того, работает ли «лучший» случайный лес с использованием «новых» предикторов лучше
model = RandomForestRegressor()
param_search = {
'n_estimators': [20, 50, 100],
'max_features': ['auto', 'sqrt', 'log2'],
'max_depth' : [i for i in range(5,15)]
}
tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)
gsearch.fit(X_train_2o, y_train_2o)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_
y_true = y_test.values
y_pred = best_model.predict(X_test)
regression_results(y_true, y_pred)
Отличные новости!! Мы значительно снизили значения RMSE и MAE, в то время как значение R-квадрат также выросло.
Feature Engineering наносит ответный удар
Давайте посмотрим, полезно ли добавление стоимости производства солнечной энергии в некотором роде для прогнозирования потребления электроэнергии.
data_consumption_2o_solar = data_consumption_2o.join(data[['Solar']]) data_consumption_2o_solar = data_consumption_2o_solar.dropna()
Сброс Train / Test + GridSearch + Проверка работоспособности
X_train_2o_solar = data_consumption_2o_solar[:'2016'].drop(['Consumption'], axis = 1)
y_train_2o_solar = data_consumption_2o_solar.loc[:'2016', 'Consumption']
X_test = data_consumption_2o_solar['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption_2o_solar.loc['2017', 'Consumption']
model = RandomForestRegressor()
param_search = {
'n_estimators': [20, 50, 100],
'max_features': ['auto', 'sqrt', 'log2'],
'max_depth' : [i for i in range(5,15)]
}
tscv = TimeSeriesSplit(n_splits=5)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)
gsearch.fit(X_train_2o_solar, y_train_2o_solar)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_
y_true = y_test.values
y_pred = best_model.predict(X_test)
regression_results(y_true, y_pred)
Вуаля, производительность модели теперь даже лучше.
График переменной важности
imp = best_model.feature_importances_
features = X_train_2o_solar.columns
indices = np.argsort(imp)
plt.title('Feature Importances')
plt.barh(range(len(indices)), imp[indices], color='b', align='center')
plt.yticks(range(len(indices)), [features[i] for i in indices])
plt.xlabel('Relative Importance')
plt.show()
Как мы видим, количество солнечной энергии не является таким сильным предсказателем энергопотребления, как другие временные предсказатели.
Финал Feature Engineering
Если вы следили за повествованием в Части 1 моего предыдущего сообщения в блоге, вы бы вспомнили, что в нашем наборе данных есть некоторый сезонный элемент, а точнее еженедельная сезонность. Таким образом, было бы разумнее использовать в качестве входных данных для модели значение потребления за неделю до указанной даты. Это означает, что если модель пытается спрогнозировать значение потребления на 8 января, ей необходимо предоставить информацию о потреблении на 1 января.
data_consumption_2o_solar_weeklyShift = data_consumption_2o_solar.copy() data_consumption_2o_solar_weeklyShift['Last_Week'] = data_consumption_2o_solar['Consumption'].shift(7) data_consumption_2o_solar_weeklyShift = data_consumption_2o_solar_weeklyShift.dropna()
Сброс Train / Test + GridSearch + Проверка работоспособности
X_train_2o_solar_weeklyShift = data_consumption_2o_solar_weeklyShift[:'2016'].drop(['Consumption'], axis = 1)
y_train_2o_solar_weeklyShift = data_consumption_2o_solar_weeklyShift.loc[:'2016', 'Consumption']
X_test = data_consumption_2o_solar_weeklyShift['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption_2o_solar_weeklyShift.loc['2017', 'Consumption']
model = RandomForestRegressor()
param_search = {
'n_estimators': [20, 50, 100],
'max_features': ['auto', 'sqrt', 'log2'],
'max_depth' : [i for i in range(5,15)]
}
tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)
gsearch.fit(X_train_2o_solar_weeklyShift, y_train_2o_solar_weeklyShift)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_
y_true = y_test.values
y_pred = best_model.predict(X_test)
regression_results(y_true, y_pred)
И мы сделали это снова .. Ошибки были уменьшены, а r-квадрат увеличился. Мы могли бы продолжать добавлять более актуальные функции, но я думаю, вы уже поняли идею!
График важности функции
Как мы правильно предположили, значение (t-7)-го дня является гораздо более сильным предиктором, чем значение (t-1)-го дня.
Заключение
В этой статье мы узнали, как моделировать данные временных рядов, проводить перекрестную проверку данных временных рядов и точно настраивать гиперпараметры нашей модели. Нам также удалось снизить RMSE с 85,61 до 54,57 для прогнозирования энергопотребления.
В Части 3 этой серии мы будем работать над тематическим исследованием, анализирующим данные временных рядов, генерируемых центрами обработки вызовов, в основном работая над анализом (пугающего) увеличения показателей отказа. Будьте на связи…
До скорого :)
