Код, стоящий за кластеризацией данных временных рядов.

В предыдущем посте мы объяснили, как работает кластеризация данных временных рядов. В этом посте мы глубоко погрузимся в сам код.

Все будет написано на python, но у большинства библиотек есть R-версия. Мы постараемся оставаться на относительно высоком уровне, но код будет содержать некоторые полезные ресурсы, если вы хотите большего.

Без лишних слов, давайте перейдем к делу.

0 - Создание данных

Чтобы облегчить интерпретацию, давайте рассмотрим теоретический пример: мы пытаемся спрогнозировать цены на золото, используя информацию с местных рынков по всему миру.

Сначала мы разрабатываем фрейм данных размером 1000 x 10. Каждая строка соответствует уникальному моменту времени, в нашем случае дню, а каждый столбец соответствует разному рынку золота. Все значения в нашем фрейме данных являются ценами.

# create data
rng = pd.date_range('2000-01-01', freq='d', periods=n_rows)
df = pd.DataFrame(np.random.rand(n_rows, n_cols), index=rng)

Приведенный выше код значительно упрощает наш пример. Одна концептуальная проблема заключается в том, что цены всегда принимают значение от 0 до 1, однако уроки кода все еще применимы.

Создав синтетический фрейм данных, давайте сделаем его грязным.

# "unclean" data
df = df.apply(lambda x: make_outliers_on_col(x), axis='index')
df = df.apply(lambda x: make_nan_on_col(x), axis='index')

Вышеупомянутые функции случайным образом вводят 10 выбросов и 10 нулевых значений в наш фрейм данных. Наш итоговый фрейм данных выглядит примерно так ...

1 - Очистка данных

Есть два основных этапа очистки данных: вменение недостающих данных и удаление выбросов. К счастью, у pandas есть несколько простых встроенных методов, которые могут нам помочь.

# interpolate missing values
df = df.interpolate(method='spline', order=1, limit=10, limit_direction='both')
# interpolate outliers
df = df.apply(lambda x: nullify_outliers(x), axis='index')
df = df.interpolate(method='spline', order=1, limit=10, limit_direction='both')

Наша стратегия здесь очень проста. Сначала мы вменяем все недостающие данные, используя сплайн-интерполяцию. Затем мы заменяем все выбросы нулевыми значениями и снова используем сплайн-интерполяцию.

В документе предложено множество методов вменения пропущенных значений, некоторые из которых включают интерполяцию (показанную выше), вычисление разложения по сингулярным значениям (SVD) и вменение K-ближайших соседей (KNN).

Если вам важна скорость, лучше всего подойдут SVD или интерполяция. KNN может обеспечить лучшие результаты, но требует больших вычислительных ресурсов.

В конце этого шага у нас будет фрейм данных, который выглядит как на рисунке 2:

2 - Кластеризация

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

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

# TS-specific features
autocorrelation = df.apply(lambda x: acf(x, nlags=3), axis='index')
partial_autocorrelation = df.apply(lambda x: pacf(x, nlags=3), axis='index')
# Signal-processing-specific features
fast_fourier_transform = df.apply(lambda x: np.fft.fft(x), axis='index')
variance = df.apply(lambda x: np.var(x), axis='index')

Отсюда мы можем выполнить кластеризацию k-средних для каждой пары наборов функций. Мы ограничиваемся двумя функциями для простоты, однако в документе упоминаются четыре потенциальных функции для обеих групп.

import numpy as np
from scipy.cluster.vq import kmeans2
# setup ts and signal features for clustering
features = [np.array([autocorrelation, partial_autocorrelation]),
            np.array([fast_fourier_transform, variance])]
for f in features:
    # cluster
    out = kmeans2(f, 2)
    cluster_centers, labels = out
    # ...

Приведенный выше код группирует каждый из наших 10 рынков золота в две отдельные группы, как показано ниже на рисунке 3.

Теперь, когда мы сгруппировали наши данные в «похожие» временные ряды », мы готовы моделировать каждую группу.

3 - Модель прогнозирования

В документе упоминается, что двунаправленный LSTM имеет наилучшую точность. Несмотря на пугающее название, двунаправленный LSTM - это всего лишь два LSTM. Первый тренируется вперед и назад с регулярными вводами. Второй обучается в обратном направлении вперед с перевернутыми входными векторами.

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

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

Прежде чем сделать это, важно отметить, что мы усреднили значения наших временных рядов в каждом кластере. Некоторые модели, такие как DeepAR, подходят для нескольких временных рядов и выводят один прогноз. Однако обычные LSTM требуют одномерных данных.

from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
# fit basic LSTM
model = Sequential()
model.add(LSTM(4, input_shape=(1, look_back)))
model.add(Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')
model.fit(trainX, trainY, epochs=n_epoch, batch_size=1, verbose=2)

Приведенный выше код итеративно запускался для всех 4 наборов данных, и их точность показана ниже на рисунке 4.

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

Рисунок 5 добавляет к нашей оценке немного большей детализации. Мы можем видеть, что функции TS работают лучше, чем функции сигналов в кластере 1, но хуже в кластере 2. В целом, средняя RMSE между кластерами 1 и 2 для каждой группы аналогична.

Теперь вам может быть интересно, почему показатели в группах так похожи. Что ж, если вы помните, наш механизм генерации данных был одинаковым для всех временных рядов ».

Цель кластеризации - выявить систематические различия в наших моделях временных рядов. Затем мы можем разработать специальную модель для каждого.

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

4 - Следующие шаги

Полный код для вышеприведенного пошагового руководства можно увидеть здесь. Однако для проекта моделирования с данными из реального мира рекомендуется попробовать больше.

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

  • Кластер по дополнительным функциям
  • Кластеризация функций TS и сигналов одновременно
  • Используйте более сложные структуры глубокого обучения
  • Ввести статические функции (в документе обсуждались архитектуры для этого)

Спасибо за внимание! Я напишу еще 30 постов, посвященных академическим исследованиям в индустрии DS. В моем комментарии есть ссылки на основной источник этого сообщения и некоторые полезные ресурсы.