После перерыва продолжаю цикл статей про одно из самых интересных направлений в статистике и науке о данных — прогнозировании временных рядов (или рядов динамики, как их первоначально называли в учебниках по эконометрике). Эта работа будет не в стиле перевода с моими комментариями, а полноценное исследование на тему эффективности прогнозных моделей: мы с вами разработаем и сравним две модели прогнозирования временных рядов — традиционную статистическую модель — реализацию модели ARIMA с сезонной компонентой и экзогенными переменными под названием SARIMAX и рекуррентную модель глубокого обучения на основе слоя LSTM. Выясним, какая их них наиболее эффективно справится с климатическими данными, которые подготовил для своей книги Франсуа Шолле «Глубокое обучение с Keras», второе издание которой вышло в 2023 году. Второе издание значительно переработано в ногу со временем, и я настоятельно рекомендую изучить эту книгу как начинающим аналитикам данных, так и бывалым представителям науки о данных с багажом знаний о временных рядах.
Попутно отвечу на накопившиеся вопросы от участников сообщества, связанных как с подготовкой данных для рекуррентных нейронных сетей, так и с объяснением деталей дальнейшего использования обученных моделей.
Приводимый код в статье обогащён моими знаниями и опробован на деле — активно пользуюсь им в проектах, связанных с применением машинного обучения, и делюсь им с вами. Но перед этим я рекомендую освежить свои знания в вопросе о том, что такое одномерные и многомерные временные ряды, а также о точечном (одношаговом) и интервальном (многошаговом) прогнозировании и их выполнении (ссылка на статью).
В конце мы рассмотрим подход в прогнозировании температуры с помощью реализации модели на основе слоя LSTM, который описывает в своей книге Франсуа, и попытаемся выяснить что в нём не так: поставленные автором цели прогнозирования не достигаются должным образом.
В общем, всем пристегнуть ремни: мы начинаем.
Для удобства повторения примеров в конце статьи будет представлена ссылка на мою записную книжку сайта kaggle.com. В её основе код с минимальными пояснениями, без объяснений, которые я привожу в данной работе. При этом я буду обновлять её по мере публикации частей исследования.
Весь представленный код в этой статье разработан в IDE Spyder и изложен в функциональном стиле. Опираясь на него, вы в большей степени расщёлкаете ту или иную задачу, связанную с таким интересным направлением в эконометрике, как прогнозирование временных рядов.
Параметры среды Python:
Python 3.10.9
matplotlib 3.8.4
numpy 1.26.4
pandas 2.2.2
statsmodels 0.14.2
keras/tensorflow 2.11.0Чтобы сократить время вычислений мы будем использовать уменьшенный набор данных, поэтому большинство примеров в данном исследовании можно выполнить с помощью ЦП. Однако тем, кто захочет провести дополнительные исследования на оригинальном наборе климатических данных я рекомендую воспользоваться google.colab.
Исследование состоит из следующих глав:
Описание набора климатических данных и его предварительная подготовка.
Разработка моделей глубокого обучения Keras с использованием слоя LSTM для аналогичной задачи.
Сравнение и оценка эффективности результатов прогнозирования.
Рассмотрение подхода прогнозирования смещённой температуры с помощью модели глубокого обучения Keras с использованием слоёв LSTM.
Вторая глава из-за большого объёма информации будет разделена на две части.
Во всех примерах статьи используются данные временных рядов о погоде, записанные на гидрометеорологической станции в Институте биогеохимии им. М. Планка (ссылка на сайт).
В этот набор данных включены замеры 14 различных метеорологических показателей (таких, как температура воздуха, атмосферное давление, влажность), выполняющиеся каждые 10 минут начиная с 2003 года. Для экономии времени мы будем анализировать данные, охватывающие период с 2009 по 2016 год, которые подготовил Франсуа Шолле для своей книги «Глубокое обучение на Python» (Deep Learning with Python by François Chollet) в обеих редакциях.
Загрузить архив с данными и распаковать его в текущий рабочий каталог можно несколькими способами. Либо с помощью утилиты get_file библиотеки tensorflow, как в книге первой редакции:
import tensorflow as tf
import os
zip_path = tf.keras.utils.get_file(
origin='https://storage.googleapis.com/tensorflow/tf-keras-datasets/jena_climate_2009_2016.csv.zip',
fname='jena_climate_2009_2016.csv.zip',
extract=True)
csv_path, _ = os.path.splitext(zip_path)Либо с помощью команд из книги Франсуа второго издания, которыми мы и воспользуемся:
!wget https://s3.amazonaws.com/keras-datasets/jena_climate_2009_2016.csv.zip
!unzip jena_climate_2009_2016.csv.zipЭтот набор данных хорошо подходит для отработки навыков и методов по прогнозированию временных рядов, поэтому советую сохранить его на локальный компьютер для быстрого доступа к нему в будущем.
Сначала импортируем необходимые модули и атрибуты, которые мы будем использовать в первой главе исследования:
from matplotlib import pyplot
from numpy import nan, where
from pandas import read_csv, DateOffsetТакже для лучшего представления графиков изменим следующие параметры метода pyplot библиотеки matplotlib на пользовательские. Размеры графиков я установил (14, 7) исходя из разрешения ноутбука, который использовал для проведения данного исследования.
pyplot.rcParams["figure.autolayout"]=True
pyplot.rcParams["figure.figsize"]=(14, 7)После того, как мы сохранили и распаковали на локальный компьютер архив с набором данных, самое время его импортировать для дальнейшего анализа. Код ниже сначала импортирует файл с помощью функции read_csv библиотеки pandas, которая в качестве первого аргумента принимает название файла с его расширением; в качестве индекса (набора идентифицирующих меток, прикреплённых к структуре данных pandas) будем использовать столбец Date Time, и чтобы pandas не импортировал значения индекса как строки, укажем ему обработать эти значения как метки даты/времени с помощью аргумента date_format.
Обращаю внимание, что в отличие от примеров прогнозирования температуры из книги Франсуа, мы будем использовать данные не в десятиминутном формате, а почасовом, поскольку это, во-первых, интереснее, чем копирование примера анализа из книги, а во-вторых, хорошенько сэкономит время выполнения дальнейших расчётов. Для этого воспользуемся методом resample() с параметром h (или «hour»): фактически выполнится группировка данных по заданному периоду времени (в нашем случае период принимается равным одному часу, который охватит шесть десятиминутных строк — от 00:00 до 00:50 включительно) и благодаря методу mean() мы вернём среднее значение для каждой группы (каждого часа).
def open_csv_and_data_resample(filename):
"""Открытие csv-файла и предварительная подготовка набора данных"""
dataframe = read_csv(filename, index_col='Date Time', date_format='%d.%m.%Y %H:%M:%S')
dataframe = dataframe.resample('h').mean()
return dataframe
df = open_csv_and_data_resample(filename='jena_climate_2009_2016.csv')Таким образом, мы импортировали csv-файл с заданными параметрами и преобразовали исходный формат набора данных с десятиминутных замеров на часовые, рассчитав среднее значение всех числовых столбцов.
Полученный набор данных имеет следующую форму:
print(df.shape)
(70129, 14)Давайте выведем на экран названия временных рядов:
print(* df.columns, sep='\n')
p (mbar)
T (degC)
Tpot (K)
Tdew (degC)
rh (%)
VPmax (mbar)
VPact (mbar)
VPdef (mbar)
sh (g/kg)
H2OC (mmol/mol)
rho (g/m**3)
wv (m/s)
max. wv (m/s)
wd (deg)Ниже приведены соответствующие им пояснения:
«Давление»,
«Температура»,
«Температура в Кельвинах»,
«Температура (точка росы)»,
«Относительная влажность»,
«Давление насыщенных паров»,
«Давление пара»,
«Дефицит давления пара»,
«Удельная влажность»,
«Концентрация водяного пара»,
«Герметичность»,
«Скорость ветра»,
«Максимальная скорость ветра»,
«Направление ветра в градусах»Взглянем на первые и последние пять строк набора данных:
print(df.head())
p (mbar) T (degC) ... max. wv (m/s) wd (deg)
Date Time ...
2009-01-01 00:00:00 996.528000 -8.304000 ... 1.002000 174.460000
2009-01-01 01:00:00 996.525000 -8.065000 ... 0.711667 172.416667
2009-01-01 02:00:00 996.745000 -8.763333 ... 0.606667 196.816667
2009-01-01 03:00:00 996.986667 -8.896667 ... 0.606667 157.083333
2009-01-01 04:00:00 997.158333 -9.348333 ... 0.670000 150.093333
[5 rows x 14 columns]print(df.tail())
p (mbar) T (degC) ... max. wv (m/s) wd (deg)
Date Time ...
2016-12-31 20:00:00 1001.410000 -2.503333 ... 1.526667 203.533333
2016-12-31 21:00:00 1001.063333 -2.653333 ... 1.250000 98.366667
2016-12-31 22:00:00 1000.511667 -3.553333 ... 1.410000 167.958333
2016-12-31 23:00:00 999.991667 -3.746667 ... 1.650000 223.600000
2017-01-01 00:00:00 999.820000 -4.820000 ... 1.960000 184.900000
[5 rows x 14 columns]По вышеприведённым таблицам можно убедиться, что период записи наблюдения составляет один час. Таким образом, в течение одного дня у нас будет 24 наблюдения, а за один обычный год (не високосный) должно накопится 8760 (24×365) наблюдения.
По конечным строкам мы видим, что к набору прицепился замер на 2017 год, от которого предлагаю избавиться, чтобы на последующих графиках данные отображались корректно. Сделаем это с помощью фильтрации pandas:
df = df[df.index.year < 2017]Выведем информацию о наборе данных вместе с метаданными:
print(df.info())
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 70128 entries, 2009-01-01 00:00:00 to 2016-12-31 23:00:00
Freq: h
Data columns (total 14 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 p (mbar) 70040 non-null float64
1 T (degC) 70040 non-null float64
2 Tpot (K) 70040 non-null float64
3 Tdew (degC) 70040 non-null float64
4 rh (%) 70040 non-null float64
5 VPmax (mbar) 70040 non-null float64
6 VPact (mbar) 70040 non-null float64
7 VPdef (mbar) 70040 non-null float64
8 sh (g/kg) 70040 non-null float64
9 H2OC (mmol/mol) 70040 non-null float64
10 rho (g/m**3) 70040 non-null float64
11 wv (m/s) 70040 non-null float64
12 max. wv (m/s) 70040 non-null float64
13 wd (deg) 70040 non-null float64
dtypes: float64(14)
memory usage: 8.0 MB
NoneНа текущий момент количество наблюдений каждого столбца набора данных составляет 70128, из них ненулевых — 70 040; разница между ними говорит о том, что среди данных есть пропущенные (пустые) значения. Также обратите внимание на частоту индексов (freq), которая определена как часовая (h).
Давайте отдельно подсчитаем для каждого столбца сумму пропущенных значений и повторяющихся строк (дубликатов):
def detect_mis_and_dub_values(dataframe, print_dubs=True):
"""Поиск пропущенных и повторяющихся значений"""
msg = "Количество пропущенных значений:"
print('\n' + msg)
print('-' * len(msg))
print(dataframe.isna().sum().sort_values(ascending=False))
if print_dubs:
msg = "Количество дубликатов:"
print('\n' + msg)
print('-' * len(msg))
print(dataframe.duplicated().sum())
detect_mis_and_dub_values(dataframe=df)Количество пропущенных значений:
--------------------------------
p (mbar) 88
T (degC) 88
Tpot (K) 88
Tdew (degC) 88
rh (%) 88
VPmax (mbar) 88
VPact (mbar) 88
VPdef (mbar) 88
sh (g/kg) 88
H2OC (mmol/mol) 88
rho (g/m**3) 88
wv (m/s) 88
max. wv (m/s) 88
wd (deg) 88
dtype: int64
Количество дубликатов:
----------------------
87Итак, у каждого временного ряда имеется 88 пустых значений. Появление дубликатов в количестве 87 штук может указывать на их связь с теми строками, на которых присутствуют одни пропущенные значения. Давайте найдём все строки с отсутствующими значениями и проверим это: сначала создадим новую переменную nan_rows с типом dataframe и с помощью фильтрации pandas присвоим ей все строки со значениями nan, а затем выведем первые и последние пять строк.
nan_rows = df[df.isnull().any(axis=1)]
print(nan_rows.head())
print(nan_rows.tail())
p (mbar) T (degC) ... max. wv (m/s) wd (deg)
Date Time ...
2014-09-24 18:00:00 NaN NaN ... NaN NaN
2014-09-24 19:00:00 NaN NaN ... NaN NaN
2014-09-24 20:00:00 NaN NaN ... NaN NaN
2014-09-24 21:00:00 NaN NaN ... NaN NaN
2014-09-24 22:00:00 NaN NaN ... NaN NaN
[5 rows x 14 columns]
p (mbar) T (degC) ... max. wv (m/s) wd (deg)
Date Time ...
2016-10-28 07:00:00 NaN NaN ... NaN NaN
2016-10-28 08:00:00 NaN NaN ... NaN NaN
2016-10-28 09:00:00 NaN NaN ... NaN NaN
2016-10-28 10:00:00 NaN NaN ... NaN NaN
2016-10-28 11:00:00 NaN NaN ... NaN NaN
[5 rows x 14 columns]Наше предположение о связи между строками с отсутствующими данными и дубликатами подтвердилось. Заметим, что пустые значения имеются в наблюдениях, зафиксированных в 2014 и 2016 годах; в предшествующих им годам (2013 и 2015 соответственно) пропусков не обнаружено.
Перед тем как обработать пропущенные значения, взглянем на распределение каждого временного ряда во времени (обращаю внимание, что отрисовка python’ом данного графика будет длительной, поскольку мы изменили настройки по умолчанию на пользовательские для лучшего представления информации):
df.plot(subplots=True, layout=(7,2))
pyplot.show()
Посмотрев на график, можно заметить, что у большинства имеющихся временных рядов наблюдаются повторяющиеся годовые колебания, что, как правило, связано с сезонными изменениями.
Итак, поскольку набор данных имеет строго расставленные метки даты и времени с частотой в один час и в данных наблюдается годовая периодичность, то предлагаю обработать отсутствующие значения по следующей логике. Обыкновенное удаление строк данных собьёт частоту и впоследствии вызовет ошибки при обучении статистических моделей прогнозирования: statsmodels будет методично выдавать предупреждения о нарушенном порядке в данных. В этой связи с помощью следующей функции мы заменим все значения nan для каждого столбца на значения, которые были зарегистрированы годом ранее.
Как и в примере с поиском пропущенных значений, мы с помощью фильтрации pandas создадим новую переменную с типом dataframe, содержащую строки со значениями nan (nan_rows). Далее с помощью цикла for пробежимся по каждому индексу (читай — временной метке) в переменной nan_rows и с помощью метода pandas DateOffset(years=1) определим для каждого найденного индекса nan_index временную метку годом ранее (prev_year_timestamp). После этого нам остаётся только заменить значения nan во всём наборе данных (либо только тех столбцах, которые будут явно указаны в аргументе функции colls) на соответствующие значения из того же времени, но годом ранее.
def replace_nan(dataframe, colls=None):
"""Поиск и замена пропущенных значений на значения предыдущего года"""
nan_rows = dataframe[dataframe.isnull().any(axis=1)]
for nan_index in nan_rows.index:
prev_year_timestamp = nan_index - DateOffset(years=1)
if colls:
dataframe.loc[nan_index, colls] = dataframe.loc[prev_year_timestamp, colls]
else:
dataframe.loc[nan_index] = dataframe.loc[prev_year_timestamp]
return dataframe
df = replace_nan(dataframe=df)Выполним команду и заново проверим набор данных на отсутствующие значения:
detect_mis_and_dub_values(dataframe=df, print_dubs=False)
Количество пропущенных значений:
--------------------------------
p (mbar) 0
T (degC) 0
Tpot (K) 0
Tdew (degC) 0
rh (%) 0
VPmax (mbar) 0
VPact (mbar) 0
VPdef (mbar) 0
sh (g/kg) 0
H2OC (mmol/mol) 0
rho (g/m**3) 0
wv (m/s) 0
max. wv (m/s) 0
wd (deg) 0
dtype: int64Пропущенные значения отсутствуют. Двигаемся дальше.
Создадим отдельную переменную с описательной статистикой и взглянем на неё:
description = df.describe()
Как по описательной статистике, так и вышеприведённому графику с распределением временных рядов во времени можно заметить, что у двух из них — «Скорость ветра» и «Максимальная скорость ветра» — среди наблюдений присутствуют выбросы в виде чудовищно отрицательных значений. Давайте построим для этих временных рядов отдельный график и взглянем на них поближе:
def show_plot_wind_colls():
"""График с 'wv (m/s)' и 'max. wv (m/s)'"""
fig, axs = pyplot.subplots(1, 2)
axs[0].plot(df['wv (m/s)'], color='orange')
axs[0].set_title('Скорость ветра')
axs[0].set_ylabel('м/с')
axs[0].set_xlabel('Дата наблюдения')
axs[1].plot(df['max. wv (m/s)'], color='green')
axs[1].set_title('Максимальная скорость ветра')
axs[1].set_xlabel('Дата наблюдения')
pyplot.suptitle('\"Скорость ветра\" и \"Максимальная скорость ветра\"')
pyplot.show()
show_plot_wind_colls()
В виду того, что скорость ветра не может быть отрицательной (по крайней мере в моём представлении) и, как уже было сказано ранее, набор данных имеет строго расставленные меткидаты и времени, то обработаем эти значения по той же логике, которую мы использовали для отсутствующих значений. Для этого мы заменим все отрицательные значения у временных рядов «Скорость ветра» и «Максимальная скорость ветра» на значения nan и впоследствии присвоим им значения, которые были зарегистрированы годом ранее.
wind_colls = ['wv (m/s)', 'max. wv (m/s)']
df[wind_colls] = df[wind_colls].apply(lambda x: where(x < 0, nan, x))Давайте посмотрим сколько пустых значений у нас после этого получилось:
detect_mis_and_dub_values(dataframe=df, print_dubs=False)
Количество пропущенных значений:
--------------------------------
wv (m/s) 4
max. wv (m/s) 4
p (mbar) 0
T (degC) 0
Tpot (K) 0
Tdew (degC) 0
rh (%) 0
VPmax (mbar) 0
VPact (mbar) 0
VPdef (mbar) 0
sh (g/kg) 0
H2OC (mmol/mol) 0
rho (g/m**3) 0
wd (deg) 0
dtype: int64У временных рядов «Скорость ветра» «Максимальная скорость ветра» наблюдается по четыре пропущенных значения. Выполним их замену на значения предыдущего года:
df = replace_nan(dataframe=df, colls=wind_colls)Таким образом, мы не только обработали все пропущенные значения, заменив их на значения предыдущего года для соответствующего временного ряда, но и разобрались с отрицательными значениями для скоростей ветра. Однако обнаруженные ранее повторяющиеся строки никуда не делись: они лишь обновили свои значения, заменив nan на замещённые значения годом ранее. Предлагаю проверить это - и вывести первые и последние пять строк переменной dubles, которую создадим с помощью уже полюбившейся нам фильтрации pandas:
dubles = df[df.duplicated()]
print(dubles.head())
print(dubles.tail())
p (mbar) T (degC) ... max. wv (m/s) wd (deg)
Date Time ...
2014-09-24 18:00:00 988.471667 15.901667 ... 1.316667 252.283333
2014-09-24 19:00:00 988.553333 15.338333 ... 1.783333 246.783333
2014-09-24 20:00:00 988.718333 15.023333 ... 1.850000 252.083333
2014-09-24 21:00:00 988.735000 14.650000 ... 2.003333 248.000000
2014-09-24 22:00:00 988.798333 13.810000 ... 1.550000 120.810000
[5 rows x 14 columns]
p (mbar) T (degC) ... max. wv (m/s) wd (deg)
Date Time ...
2016-10-28 07:00:00 988.013333 0.986667 ... 1.433333 194.750000
2016-10-28 08:00:00 988.035000 1.721667 ... 1.156667 203.750000
2016-10-28 09:00:00 988.080000 3.701667 ... 1.686667 220.550000
2016-10-28 10:00:00 987.733333 5.751667 ... 1.656667 204.766667
2016-10-28 11:00:00 987.243333 8.681667 ... 2.713333 160.683333
[5 rows x 14 columns]Что и следовало ожидать. Разумеется, мы не будем избавляться от этих повторяющихся строк, а оставим их в покое.
С учётом вышеописанной подготовки данных давайте заново построим график с распределением временных рядов во времени.
df.plot(subplots=True, layout=(7,2))
pyplot.show()
Уже другое дело!
Ещё раз отметим, что у большинства имеющихся временных рядов наблюдается выраженная годовая периодичность, с которой тесно связано такое понятие, как сезонность, анализом которой мы займёмся во второй главе исследования.
Обновим описательную статистику:
description = df.describe()
Выведем итоговую форму набора данных и информацию по нему:
print(df.shape)
(70128, 14)
print(df.info())
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 70128 entries, 2009-01-01 00:00:00 to 2016-12-31 23:00:00
Freq: h
Data columns (total 14 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 p (mbar) 70128 non-null float64
1 T (degC) 70128 non-null float64
2 Tpot (K) 70128 non-null float64
3 Tdew (degC) 70128 non-null float64
4 rh (%) 70128 non-null float64
5 VPmax (mbar) 70128 non-null float64
6 VPact (mbar) 70128 non-null float64
7 VPdef (mbar) 70128 non-null float64
8 sh (g/kg) 70128 non-null float64
9 H2OC (mmol/mol) 70128 non-null float64
10 rho (g/m**3) 70128 non-null float64
11 wv (m/s) 70128 non-null float64
12 max. wv (m/s) 70128 non-null float64
13 wd (deg) 70128 non-null float64
dtypes: float64(14)
memory usage: 10.0 MB
NoneТаким образом, набор данных представлен в виде 14-ти временных рядов, насчитывающих 70128 почасовых наблюдений, регистрация которых охватывает период с 2009-01-01 00:00:00 по 2016-12-31 23:00:00 включительно.
Как вы помните, в начале было написано про обычный год, который с учётом частоты данных в один час будет иметь общее количество наблюдений 8760 (24*365). Однако в нашем наборе климатических данных присутствуют два високосных года — 2012 и 2016, которые с учётом дополнительного дня в году будут иметь по 8784 (24*366) наблюдения. Наличие високосных годов объясняет получившуюся разницу в 48 наблюдений: 24*365*8=70080 (70080+2*24=70128).
Определение периодичности в данных — это важный шаг, поскольку он напрямую выражается в параметре сезонности m для таких статистических моделей с сезонной составляющей, как SARIMA. Изучение влияния високосных годов мы оставим за рамками данного исследования; по большому счёту, для них нужно создать дополнительные/отдельные модели. В тоже время мы не станем удалять високосные годы или их дополнительные дни из набора данных, поскольку в дальнейшем мы попробуем запустить расчёты с учётом сезонного параметра m=8760. Это должно привести к игнорированию дополнительных часов в високосные годы: модель будет учитывать только общие сезонные шаблоны для всех годов.
M=8760 … Это не ошибка! Что может быть проще для вычислений, чем сезонный период в 8760 наблюдений!? Это отличный способ проверить, сколько ваш компьютер способен продержаться, прежде чем будет уложен на лопатки математической сложностью и алгоритмической реализацией такой статистической модели, как SARIMA. Шутки шутками, но данное обстоятельство — это один из основных минусов статистических моделей с сезонной составляющей применительно к »большим» данным по причине отсутствия оптимизации расчётов. Впрочем, о сезонности, статистических методах её определения и попытках её расчёта мы более подробно поговорим во второй главе нашего исследования.
Сохраним подготовленную версию набора данных вместе с метками даты и времени в отдельный файл с расширением parquet и названием "jena_climate_2009_2016_hour_grouped_data":
df.to_parquet('jena_climate_2009_2016_hour_grouped_data.parquet')Формат parquet предпочтительнее, чем csv, поскольку по умолчанию сохраняет все метаданные набора данных (индексы, типы столбцов и т. д.) и более эффективен в плане производительности и хранения информации. К тому же, при сохранении файла в формате csv большая часть метаданных теряется (в частности частота данных — freq), поэтому при его последующем импорте необходимо явно указывать столбец с индексом и его тип/формат. Формат parquet, повторюсь, это делает по умолчанию.
А теперь отдельно посмотрим на целевой временной ряд «Температура», значения которого мы будем прогнозировать в последующих примерах исследования:
temperature = df['T (degC)']
def show_plot_temperature():
"""График температуры и её распределения"""
fig, axs = pyplot.subplots(2, 1)
axs[0].plot(temperature)
axs[0].set_title('Температура')
axs[0].set_xlabel('Дата наблюдения')
axs[0].set_ylabel('Градусов по Цельсию')
axs[0].grid()
axs[1].hist(temperature.to_numpy(), bins=50, edgecolor='black')
axs[1].set_title('Гистограмма распределения температуры')
axs[1].set_xlabel('Градусов по Цельсию')
axs[1].set_ylabel('Частота')
pyplot.show()
show_plot_temperature()
Отметим, что у температуры заметна отчётливая ежегодная периодичность, а основные значения расположены в диапазоне от -5 до 25.
Глядя на гистограмму, можно сделать предположение о нормальности распределения данных, но это не так. В рамках этого исследования мы не будем проверять рассматриваемые временные ряды на нормальность распределения. Однако применительно к целевому температурному временному ряду отмечу, что по графику «Квантиль-Квантиль» (Q-Q-график) и тесту Шапиро-Уилка данные температуры распределеныне нормально. Данное обстоятельство может привести к недостоверности результатов выполнения статистических тестов и расчёта доверительного интервала. Существуют методы, благодаря которым можно преобразовать данные временного ряда к нормальному распределению, но эти проблемы заслуживают отдельного рассмотрения в отдельной работе. В конце концов, условие нормальности входных данных не является обязательным для выполнения моделей семейства ARIMA: основное внимание для их разработки уделяется такому понятию, как стационарность временного ряда (от лат. stationarius — постоянный, неподвижный; относящийся к стоянке), с которым‑то я и предлагаю разобраться в следующей главе исследования. Там же будет рассказано о разделение рассматриваемых данных на обучающую, валидационную и тестовую части и о целях прогнозирования.
Вариант кода подготовки данных без создания графиков для командной строки терминала представлен ниже:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from numpy import nan, where
from pandas import read_csv, DateOffset
def open_csv_and_data_resample(filename):
"""Открытие csv-файла и предварительная подготовка набора данных."""
dataframe = read_csv(filename, index_col='Date Time', date_format='%d.%m.%Y %H:%M:%S')
dataframe = dataframe.resample('h').mean()
return dataframe
def replace_nan(dataframe, colls=None):
"""Поиск и замена пропущенных значений на значения предыдущего года"""
nan_rows = dataframe[dataframe.isnull().any(axis=1)]
for nan_index in nan_rows.index:
prev_year_timestamp = nan_index - DateOffset(years=1)
if colls:
dataframe.loc[nan_index, colls] = dataframe.loc[prev_year_timestamp, colls]
else:
dataframe.loc[nan_index] = dataframe.loc[prev_year_timestamp]
return dataframe
def run():
print('\n'"Выполняется открытие и преобразование файла ..."'\n')
df = open_csv_and_data_resample(filename='jena_climate_2009_2016.csv')
print("Выполняется фильтрация данных ..."'\n')
df = df[df.index.year < 2017]
print("Выполняется поиск и замена пропущенных значений ..."'\n')
df = replace_nan(dataframe=df)
wind_colls = ['wv (m/s)', 'max. wv (m/s)']
df[wind_colls] = df[wind_colls].apply(lambda x: where(x < 0, nan, x))
df = replace_nan(dataframe=df, colls=wind_colls)
print("Сохранение подготовленного набора данных в отдельный файл ..."'\n')
df.to_parquet('jena_climate_2009_2016_hour_grouped_data.parquet')
print("... завершено"'\n')
if __name__ == "__main__":
run()
Dynasaur
Прогнозирование такого примера доступно и ребёнку без всяких сложных моделей. Взяв карандаш не сложно продолжить волны сколь угодно вперёд - зимой холоднее, летом жарче. И результат будет не менее точным. Можете привести пример удачного прогнозирования, где простому взгляду на график не так очевидно?