Регрессия больших данных — Полный конвейер — Часть 1/2

TL; DR

В этой статье, состоящей из двух частей, представлен всесторонний взгляд на реальную проблему больших данных регрессии, включая код и подробный анализ. Часть 1 (эта статья) в основном посвящена EDA и моделированию поверхностного обучения, а часть 2 в основном посвящена моделированию глубокого обучения.

Повестка дня

1. Описание проблемы и набора данных
2. Исследовательский анализ данных (EDA)
3. Моделирование поверхностными методами

1. Описание проблемы и набора данных

Первоначальная проблема была обобщена для соблюдения конфиденциальности клиента. 500 различных датчиков производят выборку данной сцены с частотой 1 Гц (каждую секунду) и ссылаются на 2 числовые метки (представьте себе, что это производительность, качество и т. д.). Запрошенный результат состоит в том, чтобы предсказать эти 2 метки на основе сенсорных данных.

Описание набора данных:

• var0 ➔ Отметка времени
• var1 … var501 ➔ Сенсорные образцы (условные значения с отсутствующими значениями)
• результат 502, результат 502 ➔ Ярлыки (должны быть предсказаны)
• ~3 миллиона строк, с 01.09.2021 по 04.10.2021 (~1 месяц)

Данные хранятся в огромном (~7 Гб) CSV-файле, и для их эффективной загрузки применяется Dask:

import dask.dataframe as dd
data_df = dd.read_csv(path.join(base_path, 'raw_data.csv'))
data_df = data_df.drop('Unnamed: 0', axis=1)
display(data_df.head())

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

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

!pip install sweetviz dask[dataframe] dask-ml dask_xgboost lazypredict==0.2.9
import pickle
import numpy as np
import pandas as pd
from os import path
import seaborn as sns
import sweetviz as sv
from glob import glob
import dask.dataframe as dd
import matplotlib.pyplot as plt
from sklearn.utils import shuffle
from dask.distributed import Client
from sklearn.metrics import r2_score
from sklearn.pipeline import Pipeline
from dask.diagnostics import ProgressBar
from dask_ml.xgboost import XGBRegressor
from dask_ml.decomposition import IncrementalPCA
from dask_ml.preprocessing import StandardScaler
from lazypredict.Supervised import LazyRegressor
from sklearn.model_selection import train_test_split
from dask_ml.model_selection import train_test_split as dd_train_test_split

2. Исследовательский анализ данных (EDA)

2.1. Тип данных

Изучите каждый тип данных столбца, применив метод describe dask:

display(data_df.describe())

Можно заметить, что все столбцы имеют тип float64.

2.2. Количество образцов

Исследуйте количество образцов, вычислив форму кадра данных Dask:

with ProgressBar():
    print("Data contains %d samples" % data_df.shape[0].compute())

Выполнение приведенного выше кода вернуло отчет о 2923189 выборках, что действительно соответствует примерно 1 месяцу выборки с частотой 1 Гц (~ 31x24x60x60).

2.3. Проверка этикеток

Нанесите на график две метки, обозначенные как «результат 502» и «результат 503»:

results_list = ['var0', "result 502", "result 503"]
with ProgressBar():
   results_df = data_df[results_list].compute().set_index('var0')
fig, axes = plt.subplots(2, 1, figsize=(20,12))
results2_df = results_df.iloc[:int(1e5),:]
results_df.plot(title='Labels (raw) - full', grid=True, ax=axes[0])
results2_df.plot(title='Labels (raw) - zoom', grid=True, ax=axes[1])
plt.show()

Выполнение приведенного выше кода над данными приводит к следующим цифрам:

Можно заметить, что две метки коррелируют, в то время как «результаты 503», по-видимому, действуют прерывистым образом, т. е. с периодическим падением до нуля. Это наблюдение подразумевает, что, возможно, стоит проверить, можно ли разделить данные (выборки) на 2 группы в зависимости от статуса «результаты 503» ➔ см. раздел 2.6.

2.4. Контроль образцов

Проверьте количество пропущенных значений (None, NaN, NaT):

missing_values = data_df.isnull().sum()
percent_missing = ((missing_values / data_df.index.size) * 100)
with ProgressBar():
   percent_missing_s = percent_missing.compute()
print('Missing Values mean: %.2f%%' % percent_missing_s.mean())
title_str = 'Missing Values [%]'
percent_missing_s.plot(title=title_str, grid=True, figsize=(20,5))

Выполнение приведенного выше кода над данными привело к следующему рисунку:

Отсутствующие значения означают: 15,71 %

Ниже приводится «увеличение» для нескольких репрезентативных выборок (код + цифры), чтобы лучше понять поведение пропущенных значений:

with ProgressBar():
   few_features_list = ['var100', 'var200', 'var300', 'var400', 'var500']
   few_features_df = data_df[['var0'] + few_features_list]
   few_features_df = few_features_df.compute().set_index('var0')
fig, axes = plt.subplots(5, 1, figsize=(20,25))
for k, features in enumerate(few_features_df.iteritems()):
   title_str = '%s (raw)' % features[0]
   features[1].plot(title=title_str, grid=True, ax=axes[k])
plt.show()

Можно заметить, что пропущенные значения представляют собой большие отбрасываемые сегменты в соответствующих сигналах, что означает, что следует рассмотреть стратегию заполнения нулем. Рациональность заключается в том, что «отсутствующие данные» можно рассматривать как «время простоя», и, следовательно, нулевое заполнение кажется подходящим выбором.

Обратите внимание, что иногда имеет смысл использовать другие стратегии, такие как заполнение предыдущим, средним, медианным и т. д. В некоторых случаях может иметь смысл отбрасывать признаки с отсутствующими значениями. Следующий код может быть полезен при применении такой стратегии (удаление столбцов с более чем 50% отсутствующих значений):

columns_to_drop = percent_missing_s[percent_missing_s >= 50].index
data_cleaned_df = data_df.drop(list(columns_to_drop), axis=1)
print('Columns number reduced from %d to %d after clean-up' % \
      (len(data_df.columns), len(data_cleaned_df.columns)))

2.5. Корреляционный анализ

Коэффициент корреляции Пирсона (r) измеряет статистическую взаимосвязь (или ассоциацию) между двумя непрерывными переменными, основанную на ковариации. Он находится в диапазоне от -1 до +1 и количественно определяет направление и силу линейной связи между двумя переменными.
Эванс (1996) предложил следующее для абсолютного значения:

• 0,0 ‹= r ‹ 0,2 ➔ «очень слабый»
• 0,2 ‹= r ‹ 0,4 ➔ «слабый»
• 0,4 ‹= r ‹ 0,6 ➔ «умеренная»
• 0,6 ‹= r ‹ 0,8 ➔ «сильный»
• 0,8 ‹= r ‹ 1,0 ➔ «очень сильный»

Ниже приведена Correlation-Matrix (code+figure), сгенерированная методом Dask corr и построенная с помощью метода matshow библиотеки matplotlib:

with ProgressBar():
   corr_pearson_df = data_df.corr(method='pearson').compute()
matfig = plt.figure(figsize=(10,10))
plt.matshow(corr_pearson_df.corr(), fignum=matfig.number)
plt.show()

Можно заметить, что в корреляционной матрице мало желтых пятен. Это подразумевает наличие взаимной корреляции между выборками и возможный потенциал их объединения в объединенные объекты и/или рассмотрение методов уменьшения размерности (например, PCA, t-SNE и т. д.) для использования взаимной информации между выборками ➔ см. раздел 3.1.

Далее интересно сосредоточиться на последних двух строках/столбцах, которые представляют корреляцию между образцами и метками:

CORRELATION_THR = 0.4
correlation_dict = {}
fig, axes = plt.subplots(2, 1, figsize=(20,10))
for k, label_type in enumerate(('502', '503')):
   corr_s = corr_pearson_df.loc['result %s' % label_type]
   corr_s = correlation_s.drop(['result 502', 'result 503'])
   print('\nSignificant Feature for results %s (>%.1f):' % \
         (label_type, CORRELATION_THR))
   print(corr_s[corr_s > CORRELATION_THR].index.values)
   title_str = 'Result %s - Pearson Correlation' % label_type
   corr_s.plot(ax=axes[k], grid=True, title=title_str)
   correlation_dict[label_type] = corr_s
plt.show()

Можно заметить, что обе метки имеют «умеренную» корреляцию и выше (0,4+) среди нескольких образцов, что обнадеживает.

Следующий код вычисляет и сохраняет функции с высокой степенью корреляции:

fig, axes = plt.subplots(2, 1, figsize=(20,15))
high_correlated_dict = {}
for k, label_type in enumerate(('502', '503')):
   corr_s = correlation_dict[label_type]
   columns_to_keep = ['var0', 'result 502', 'result 503'] + \
                     corr_s[corr_s > CORRELATION_THR].index.tolist()
   with ProgressBar():
      df = data_df[columns_to_keep].compute().set_index('var0')
   print('\nHigh Correlated features %s:' % str(df.shape))      
   display(df.head())
   title_str = 'High Correlated Features (result %s)' % label_type
   plot_df = corr_pearson_df.loc[df.columns]['result %s'%label_type]
   plot_df = plot_df.drop(['result 502', 'result 503'])
   plot_df.plot(kind='bar', title=title_str, grid=True, ax=axes[k])
   high_correlated_dict[label_type] = df
plt.show()

1-я метка («результат 502») оказалась с 8 высококоррелированными образцами, а 2-я метка («результат 503») — 34 высококоррелированными образцами.

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

2.6. Активный против неактивного

Этот раздел основан на предыдущих разделах. Изучение двух меток в разделе 2.3 предполагает, что данные (выборки), возможно, можно разделить на 2 группы в зависимости от статуса «results 503», т. е. «Активно» или «Неактивно». Кроме того, в разделе 2.5 подразумевается, что только небольшая часть выборок коррелирует с метками, поэтому может быть полезно сосредоточиться на этих функциях (с точки зрения вычислений и памяти, в основном для анализа больших данных).

Следующий код применяет функцию compare_intra SweetViz для проведения такого анализа. Эта функция принимает логический ряд в качестве одного из аргументов, а также явный кортеж name для именования результирующих наборов данных (true, false). Обратите внимание, что внутренне это создает 2 отдельных кадра данных для представления каждой результирующей группы, в данном случае Активной и Неактивной.

for k, label_type in enumerate(('502', '503')):
   
   html_log = 'report_%s.html' % label_type
   feature_config = sv.FeatureConfig(
      force_num=high_correlated_df.columns.tolist()
   )
   report = sv.compare_intra(high_correlated_df,
                             high_correlated_df['result 503'] == 0,
                             ['Stop', 'Active'],
                             'result %s' % label_type,
                             feature_config,
                             pairwise_analysis='auto')
   report.show_html(html_log)

Выполнение приведенного выше кода приводит к созданию двух HTML-отчетов, по одному на каждую метку («502» и «503»). Каждый отчет предоставляет подробное представление о каждом образце, а также мощное представление парного анализа (действующее как корреляционная матрица).

Приведенный выше рисунок взят из двух отчетов HTML. Левая сторона относится к отчету "результат 502", а правая сторона относится к отчету "результат 503". Верхний ряд обеспечивает высокоуровневое статистическое сравнение между двумя группами. Средняя часть посвящена одной выборке с наиболее высокой корреляцией ("var109'), а нижняя часть представляет матрицы парной корреляции для каждой группы. Обратите внимание, что SweetViz compare_intra плохо работает с большим количеством функций, поэтому важно применять его к выбранному подмножеству, например. наиболее коррелированные.

2.7. Уменьшение размерности

Таким образом, анализ основных компонентов (PCA) применяется для использования максимальной дисперсии при сохранении наименьшего количества компонентов. Dask Incremental PCA применяется для работы с большими объемами данных в условиях ограничения памяти.

lr = Pipeline(steps=[
     ('scale', StandardScaler()),
     ('pca', IncrementalPCA()),
])
dropped_list = ['var0', 'result 502', 'result 503']
data_filt_df = data_df.drop(dropped_list, axis=1).fillna(0)
X_arr = data_filt_df.to_dask_array(lengths=True).astype(float)
lr_fitted = lr.fit(X_arr)
pca = lr_fitted['pca']
print('PCA explained variance (%d components):' % pca_components_num)
print(pca.explained_variance_ratio_)
print(pca.singular_values_)
explained_variance = 0.9
cumsum_var = pca.explained_variance_ratio_.cumsum() > explained_variance
idx = cumsum_var.argmax()
print('Number of components needed for having at least %.2f is equal to %d' % (explained_variance, idx))
pd.Series(pca.explained_variance_ratio_.cumsum()).plot(title='PCA explained variance (n=%d)' % pca_components_num, grid=True, figsize=(20,5))

Можно заметить, что для объяснения 90% дисперсии требуется 121 компонент. В настоящее время мы используем PCA для всех данных для целей EDA. В дальнейшем PCA с моделированием можно применять следующим образом:

• Подбор PCA только на тренировочном наборе ➔ pca.fit(X_train)
• Примените сопоставление (преобразование) как к обучающему набору, так и к тестовому набору ➔ X_train = pca.transform(X_train), X_test = pca.transform(X_test)

3. Моделирование поверхностными методами

3.1. Особенности Инжиниринг

EDA, представленный в разделе 2, обрабатывал выборки ('var1', …, 'var501') как функции. На этом этапе при работе с Feature Engineering можно рассмотреть следующую схему (можно применить одну или несколько):

for label_type in ('502', '503'):
   # only small portion, due to memory limitations:
   df = high_correlated_dict[label_type].iloc[:int(1e4),:]
   # Fill missing values with zeros:
   df = df.fillna(0)
   # Structure the data (+shuffle):
   X, y = shuffle(df.drop(['result 502', 'result 503'], axis=1),
                  df['result %s' % label_type],
                  random_state=42)
   # Split data into train and test sets (80:20 ratio):
   cv_tuple = train_test_split(X, y, test_size=.2, random_state=0)
   X_train, X_test, y_train, y_test = cv_tuple
   # Create an instance of the estimator, and fit it to the data:
   reg = LazyRegressor(verbose=0,
                       ignore_warnings=False,
                       custom_metric=None,
                       predictions=True)
   lazy_models, lazy_predictions = reg.fit(X_train, 
                                           X_test,
                                           y_train,
                                           y_test)
   # Results analysis:
   print('Lazy Models (label = result %s)' % label_type)
   display(lazy_models)
   if not lazy_models.empty:
      fig, axs = plt.subplots(1, 1, figsize=(20, 5), sharex=True)
      sns.set_theme(style="whitegrid")
      sns.barplot(x=lazy_models.index,
                  y="R-Squared", 
                  data=lazy_models,
                 ax=axs)
      axs.set(ylim=(0, 1))
      axs.set(title='Lazy Models (label = result %s)' % label_type)
      plt.xticks(rotation=90)
      plt.show()

client = Client(processes=True, threads_per_worker=1)
#client = Client('scheduler-address:8786')
for label_type in ('502', '503'):
   # Structure the data (high-correlated features only):
   df = high_correlated_dict[label_type]
   high_correlated_dd = dd.from_pandas(df , npartitions=3)
   X = high_correlated_dd.drop(['result 502', 'result 503'], axis=1)
   y = high_correlated_dd['result %s' % label_type]
   # Split data into train and test sets (80:20 ratio):
   cv_tpl = dd_train_test_split(X, y, test_size=.2, random_state=0)
   X_train, X_test, y_train, y_test = cv_tpl
   # Fit the model
   est = XGBRegressor()
   est.fit(X_train, y_train)
   # Results analysis:
   test_df = pd.DataFrame(y_test.compute())
   test_df['predicted'] = est.predict(X_test).compute()
   accuracy_r2 = r2_score(test_df['result %s' % label_type],
                          test_df['predicted'])
   print('accuracy=%.2f' % accuracy_r2)
   feature_import_df = pd.DataFrame(index=X_test.columns,
                                    data=est.feature_importances_)
   title_str = 'Result %s - Features Import. (XGboost)' % label_type
   feature_import_df.plot(title=title_str, 
                          kind='bar',
                          legend=None,
                          figsize=(20,5))
   fig, axes = plt.subplots(2, 1, figsize=(20,12))
   title_str = 'Result %s - Act. vs. Predicted - full' % label_type
   test_df.plot(title=, grid=True, ax=axes[0])
   title_str = 'Result %s - Act. vs. Predicted - first 10K samples'    
   test_df_zoom = test_df.iloc[:int(1e4),:]
   test_df_zoom.plot(title=title_str, grid=True, ax=axes[1])
   plt.show()