Поиск аномалий временного ряда при помощи XGboost
Нужно найти аномалии на временном ряде и для этой задачи выбрал подход машинного обучения, а именно XGboost.
Есть датасет, состоящей из 2-х параметров: секунд и значения ряда(децибелы).
Проблема: y_test, получаемая из функции PrepareData() почему-то пустой и никак не пойму почему так происходит. У Вас есть идеи для решения проблемы?
Код:
def code_mean(data, cat_feature, real_feature):
"""
Возвращает словарь, где ключами являются уникальные категории признака cat_feature,
а значениями - средние по real_feature
"""
return dict(data.groupby(cat_feature)[real_feature].mean())
Для удобства делаю агрегацию
df1 = pd.read_csv("data.csv", sep=',')
df1.columns = ['Time, seconds', 'y']
df1['Time, seconds'] = df1['Time, seconds']*10000
df1['Time, seconds'] = pd.to_datetime(df1['Time, seconds'], unit='s')
df1.set_index('Time, seconds', inplace=True)
aggregated_df1 = df1.resample('H').mean()
data = pd.DataFrame(aggregated_df1)
data.index = pd.to_datetime(data.index)
data["hour"] = data.index.hour
data["weekday"] = data.index.weekday
data['is_weekend'] = data.weekday.isin([5,6])*1
data.head()
Функция, которая сразу же будет возвращать разбитые на трейн и тест датасеты и целевые переменные. Тут-то и все проблемы:
def prepareData(data, lag_start=5, lag_end=20, test_size=0.15):
data = pd.DataFrame(data.copy())
data.columns = ["y"]
# считаем индекс в датафрейме, после которого начинается тестовыый отрезок
test_index = int(len(data)*(1-test_size))
# добавляем лаги исходного ряда в качестве признаков
for i in range(lag_start, lag_end):
data["lag_{}".format(i)] = data.y.shift(i)
data.index = pd.to_datetime(data.index)
data["hour"] = data.index.hour
data["weekday"] = data.index.weekday
data['is_weekend'] = data.weekday.isin([5,6])*1
# считаем средние только по тренировочной части, чтобы избежать лика
data['weekday_average'] = data['weekday'].apply(lambda x: code_mean(data[:test_index], 'weekday', "y").get(x))
data["hour_average"] = data['hour'].apply(lambda x: code_mean(data[:test_index], 'hour', "y").get(x))
# выкидываем закодированные средними признаки
data.drop(["hour", "weekday"], axis=1, inplace=True)
data = data.dropna()
data = data.reset_index(drop=True)
# разбиваем весь датасет на тренировочную и тестовую выборку
X_train = data.loc[:test_index].drop(["y"], axis=1)
y_train = data.loc[:test_index]["y"]
X_test = data.loc[test_index:].drop(["y"], axis=1)
y_test = data.loc[test_index:]["y"]
return X_train, X_test, y_train, y_test
Ну и основная функция для прогнозирования:
def XGB_forecast(data, lag_start=5, lag_end=20, test_size=0.15, scale=1.96):
# исходные данные
X_train, X_test, y_train, y_test = prepareData(aggregated_df1, lag_start, lag_end, test_size)
dtrain = xgb.DMatrix(X_train, label=y_train)
dtest = xgb.DMatrix(X_test)
# задаём параметры
params = {
'objective': 'reg:linear',
'booster':'gblinear'
}
trees = 1000
# прогоняем на кросс-валидации с метрикой rmse
cv = xgb.cv(params, dtrain, metrics = ('rmse'), verbose_eval=False, nfold=10, show_stdv=False, num_boost_round=trees)
# обучаем xgboost с оптимальным числом деревьев, подобранным на кросс-валидации
bst = xgb.train(params, dtrain, num_boost_round=cv['test-rmse-mean'].argmin())
# можно построить кривые валидации
#cv.plot(y=['test-mae-mean', 'train-mae-mean'])
# запоминаем ошибку на кросс-валидации
deviation = cv.loc[cv['test-rmse-mean'].argmin()]["test-rmse-mean"]
# посмотрим, как модель вела себя на тренировочном отрезке ряда
prediction_train = bst.predict(dtrain)
plt.figure(figsize=(15, 5))
plt.plot(prediction_train)
plt.plot(y_train)
plt.axis('tight')
plt.grid(True)
# и на тестовом
prediction_test = bst.predict(dtest)
lower = prediction_test-scale*deviation
upper = prediction_test+scale*deviation
if len(y_test) > 0:
Anomalies = np.array([np.NaN]*len(y_test))
Anomalies[y_test<lower] = y_test[y_test<lower]
else:
Anomalies = np.array([])
plt.figure(figsize=(15, 5))
plt.plot(prediction_test, label="prediction")
plt.plot(lower, "r--", label="upper bond / lower bond")
plt.plot(upper, "r--")
plt.plot(list(y_test), label="y_test")
plt.plot(Anomalies, "ro", markersize=10)
plt.legend(loc="best")
plt.axis('tight')
plt.title("XGBoost Mean absolute error {} users".format(round(mean_absolute_error(prediction_test, y_test))))
plt.grid(True)
plt.legend()
Вызов:XGB_forecast(aggregated_df1, test_size=0.2, lag_start=5, lag_end=30)