Поиск аномалий временного ряда при помощи 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)


Ответы (0 шт):