Машинное обучение с CNTK от Microsoft: анализ временных рядов

Эта часть будет посвящена регрессионному анализу, если быть точным — анализу временных рядов. Чем эта задача принципиально отличается от задачи классификации, которую мы рассматривали в предыдущих 2 частях (первая, вторая)? Тем, что мы исследуем не принадлежность какого-либо объекта классу (собака, кошка и т.д.), а пытаемся построить функцию y = f(x), которая вычисляет значение переменной y (зачастую количественной, например, частота звука, стоимость и т.д., иными словами, — что-то, что можно измерить) в зависимости от входных параметров.

Многие из вас не раз видели подобную картинку:

Это обычная линейная регрессия, даже не нейронная сеть, хотя они имеют много общего, можно даже сказать, что это нейронная сеть всего из одного нейрона и без функции активации. Что это значит? Это значит, что, работая с нейронной сетью, сразу можно выявить нелинейные зависимости в наших данных. Одним из примеров выявления нелинейной зависимости является тот же XOR, который строили в первой части. Это «базовые» преимущества нейронной над обычной моделью регрессии.

Помимо этого рассмотрим также предварительную обработку данных при регрессионном анализе и разберемся, зачем она нужна. Для примера возьмем «шаблонный» набор данных Auto MPG, который содержит характеристики автомобилей 70–80-х годов. Исходя из них нужно оценить расход топлива автомобиля при езде в городском цикле. Проводя регрессионный анализ, очень важно понимать, какого типа у нас переменные и под какие типы распределений они подпадают. Пока что импортируем все необходимые модули, загрузим датасет и посмотрим, что там внутри:

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import cntk as C

data = pd.read_csv('Data/auto_mpg.csv')

data.hist(column=['mpg','cylinders','displacement','horsepower','weight','acceleration','model_year'])
plt.show()

Почему-то нет гистограммы для лошадиных сил (horsepower), почему же? Потому что присутствуют не все значения. Отсутствие и неполнота данных — это вполне нормальное явление, пока для наглядности заменим эти значения на 0 в csv файле.

Теперь по очереди: очевидно, что поле acceleration имеет нормальное распределение, cylinders больше похоже на дискретное, но неравномерное распределение, поэтому над ним будет сложно сделать какое-либо преобразование.

displacement похоже на логнормальное распределение с просадкой между 250 и 300. Потому что не так много двигателей в то время имело объем между 250 и 300 кубическими дюймами (4-4.9 литра).

horsepower — точно так же логнормальное распределение, отсутствующие значения (которые временно заменили 0) можно заменить средним значением, либо же найти соответствующие характеристики автомобилей исходя из объема/года выпуска и т.д.

model_year — то же самое, что и cylinders. mpg и weight, которые имеют логнормальное распределение.

К чему мы сейчас описывали данные так подробно? К тому, что мы будем обучать не одну, а две сети сразу и сравним их в действии. Одна будет учиться на «голых» (raw в англоязычной литературе) данных, вторая получит на вход преобразованные в соответствии с распределениями, описанными выше. Во многих туториалах и гайдах датасет нормализуют или стандартизируют «не глядя», тем самым «ломая» свойства и природу первоначальных распределений.

Попробуем сеть с такой конфигурацией:

def create_model(input_):
   with C.layers.default_options(activation=C.selu, init=C.glorot_uniform()):
       model = C.layers.Sequential([
           C.layers.Dense(25),
           C.layers.Dense(25),
           C.layers.Dense(OUTPUT_DIM, activation=None)
       ])

   return model(input_)

Функция активации — SELU (подробнее читаем здесь, очень хорошо описаны свойства каждой из функций, будет полезно «продвинутым» data scientist’ам). CNTK содержит большинство функций активаций из этого списка. Если вы захотите реализовать свою функцию, вам придется не просто создать Function, а работать с UserFunction, вот пример из официальной документации.

Гиперпараметры:

EPOCHS_COUNT = 2000
LEARNING_RATE = 0.02
INPUT_DIM = 7
OUTPUT_DIM = 1

Теперь давайте разобьем данные на тренировочные (70%), тестовые (15%), валидационные (15%) для обычной сети (без предварительных преобразований):

probs = np.random.rand(len(data))
training_mask = probs < 0.7
test_mask = (probs>=0.7) & (probs < 0.85)
validation_mask = probs >= 0.85

df_training = data[training_mask]
df_test = data[test_mask]
df_validation = data[validation_mask]

Этой функцией будем преобразовывать дата-фрейм в массив, который сможет понять CNTK:

def get_df_data(df):
   X = np.array(df[['cylinders', 'displacement', 'horsepower', 'weight', 'acceleration', 'model_year', 'origin']], dtype='float32')
   y = np.array(df['mpg'], dtype='float32')
   y = y.reshape(-1, 1)

   return X, y

Далее все как и в предыдущей части: конфигурируем модель, «подвязываем входы и выходы», проводим обучение. Поскольку это модель регрессионного анализа, целесообразно в качестве функции ошибки взять Squared Error. Метод оптимизации — Adam:

#«Подвязка» для входа
input_ = C.input_variable(INPUT_DIM, 'float32')
#«Подвязка» для выхода
output_ = C.input_variable(OUTPUT_DIM, 'float32')

#Нейронная сеть
z = create_model(input_)

#«Расписание» обучения
lr_schedule = C.learning_rate_schedule(LEARNING_RATE, C.UnitType.minibatch)

learner = C.adam(z.parameters, lr_schedule, momentum=0.9)
loss = C.squared_error(z, output_)
error = C.squared_error(z, output_)

trainer = C.Trainer(z, (loss, error), [learner])

train_X, train_y = get_df_data(df_training)
test_X, test_y = get_df_data(df_test)
val_X, val_y = get_df_data(df_validation)

datamap = {
   input_: train_X,
   output_: train_y
}


for i in range(EPOCHS_COUNT):
   trainer.train_minibatch(datamap)
   if i % 100 == 0:
       print('Epoch: {epoch} Loss: {val}'.format(val=trainer.previous_minibatch_loss_average, epoch=i))

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

def pre_process_data(df_train, df_test, df_val):
   acc_m = df_train['acceleration'].mean()
   acc_std = df_train['acceleration'].std()

   cyl_min = df_train['cylinders'].min()
   cyl_max = df_train['cylinders'].max()

   hp_m = df_train['horsepower'].mean()
   hp_std = df_train['horsepower'].std()

   model_year_min = df_train['model_year'].min()
   model_year_max = df_train['model_year'].max()

   df_train['acceleration'] = (df_train['acceleration'] - acc_m) / acc_std
   df_test['acceleration'] = (df_test['acceleration'] - acc_m) / acc_std
   df_val['acceleration'] = (df_val['acceleration'] - acc_m) / acc_std

   df_train['horsepower'] = (df_train['horsepower'] - hp_m) / acc_std
   df_test['horsepower'] = (df_test['horsepower'] - hp_m) / acc_std
   df_val['horsepower'] = (df_val['horsepower'] - hp_m) / acc_std

   df_train['displacement'] = np.log(df_train['displacement'])
   df_test['displacement'] = np.log(df_test['displacement'])
   df_val['displacement'] = np.log(df_val['displacement'])

   df_train['mpg'] = np.log(df_train['mpg'])
   df_test['mpg'] = np.log(df_test['mpg'])
   df_val['mpg'] = np.log(df_val['mpg'])

   df_train['weight'] = np.log(df_train['weight'])
   df_test['weight'] = np.log(df_test['weight'])
   df_val['weight'] = np.log(df_val['weight'])

   return df_train, df_test, df_val



df_training, df_test, df_validation = pre_process_data(df_training, df_test, df_validation)

processed_train_X, processed_train_y = get_df_data(df_training)
processed_test_X, processed_test_y = get_df_data(df_test)
processed_val_X, processed_val_y = get_df_data(df_validation)

Конфигурируем новую сеть:

input_ = C.input_variable(INPUT_DIM, 'float32')
output_ = C.input_variable(OUTPUT_DIM, 'float32')

z_norm = create_model(input_)

lr_schedule = C.learning_rate_schedule(LEARNING_RATE, C.UnitType.minibatch)

learner = C.adam(z_norm.parameters, lr_schedule, momentum=0.9)
loss = C.squared_error(z_norm, output_)
error = C.squared_error(z_norm, output_)

trainer = C.Trainer(z_norm, (loss, error), [learner])

Проводим обучение:

datamap = {
   input_: processed_train_X,
   output_: processed_train_y
}


for i in range(EPOCHS_COUNT):
   trainer.train_minibatch(datamap)
   # error_function.append(trainer.previous_minibatch_loss_average)
   if i % 100 == 0:
       print('Epoch: {epoch} Loss: {val}'.format(val=trainer.previous_minibatch_loss_average, epoch=i))

Теперь самый интересный момент. Сравниваем результаты сетей:

test_scores = z.eval(test_X)
val_scores = z.eval(val_X)
processed_test_scores = z_norm.eval(processed_test_X)
processed_val_scores = z_norm.eval(processed_val_X)

test_scores = test_scores.flatten()
test_y = test_y.flatten()
val_scores = val_scores.flatten()
val_y = val_y.flatten()

processed_test_scores = processed_test_scores.flatten()
processed_test_y = processed_test_y.flatten()
processed_val_scores = processed_val_scores.flatten()
processed_val_y = processed_val_y.flatten()

Для наглядности отсортируем значения по возрастанию, для того чтобы было подобие прямой линии (когда увидите график, станет понятнее, зачем это, уже совсем скоро):

_ = val_y.argsort()
val_scores = val_scores[_]
val_y = val_y[_]

_ = processed_test_y.argsort()
processed_test_scores = processed_test_scores[_]
processed_test_y = processed_test_y[_]

_ = processed_val_y.argsort()
processed_val_scores = processed_val_scores[_]
processed_val_y = processed_val_y[_]


fig = plt.figure(figsize=(7,7))
plt.subplot(2, 2, 1, aspect=1)
plt.plot(test_scores, 'y^')
plt.plot(test_y, 'bo')
plt.xlim(-5, 65)
plt.ylim(0, 50)
plt.title('Test Raw')

plt.subplot(2, 2, 2, aspect=1)
plt.plot(val_scores, 'y^')
plt.plot(val_y, 'bo')
plt.xlim(-5, 65)
plt.ylim(0, 50)
plt.title('Validation Raw')

plt.subplot(2, 2, 3, aspect=1)
# Преобразование обратное функции логарифмирования 
# (см. функцию preprocess_data)
plt.plot(np.exp(processed_test_scores), 'y^')
plt.plot(np.exp(processed_test_y), 'go')
plt.xlim(-5, 65)
plt.ylim(0, 50)
plt.title('Test Preprocessed')

plt.subplot(2, 2, 4, aspect=1)
# Преобразование обратное функции логарифмирования 
# (см. функцию preprocess_data)
plt.plot(np.exp(processed_val_scores), 'y^')
plt.plot(np.exp(processed_val_y), 'go')
plt.xlim(-5, 65)
plt.ylim(0, 50)
plt.title('Validation Preprocessed')

plt.show()

Получается следующая картина (синим цветом показаны результаты сети без предварительной обработки на тестовых и валидационных данных, зеленым — с предварительной обработкой):

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

Давайте добавим в первый слой еще 25 нейронов:

def create_model(input_):
   with C.layers.default_options(activation=C.selu, init=C.glorot_uniform()):
       model = C.layers.Sequential([
           C.layers.Dense(50),
           C.layers.Dense(25),
           C.layers.Dense(OUTPUT_DIM, activation=None)
       ])

   return model(input_)

Уже лучше, сети практически сравнялись в точности своей работы. Почему так происходит и что из этого следует? Это проблема сходимости методов градиентного спуска и обусловленность Гессиана, порождаемая нашими данными и сетью. Это весьма продвинутая в плане изучения тема, и если вы хотите понять/разобраться, что вы сейчас прочитали, то рекомендуем также прочитать вот эту книгу.

Продолжение следует…

Авторы:
Александр Ганджа, CTO DataTrading
Богдан Домненко, Data Scientist DataTrading

Ещё интересное для вас:
Тест: какой язык программирования вам стоит выбрать для изучения?
Тест: как хорошо вы разбираетесь в Data Science?
Соревнования и бесплатная онлайн-школа для программистов