Ръководство стъпка по стъпка за проектиране на функции за многомерни времеви серии

Добавяне на нови функции въз основа на обобщена статистика с помощта на Python

Това е ръководство стъпка по стъпка за проектиране на функции за „прогнозиране на многовариантни времеви редове“. Ще научите как да изчислявате няколко подвижни статистики. Добавянето им към обяснителните променливи често води до по-добро прогнозиране.

Въведение

Авторегресия

Многовариантният времеви ред съдържа две или повече променливи. Вижте по-долу за пример. Често тези набори от данни се изучават с цел прогнозиране на една или повече от тези променливи.

Повечето модели за прогнозиране се основават на авторегресия. Това се равнява на решаване на задача за регресия на контролирано обучение. „Бъдещите стойности на серията са целевите променливи. Входните обяснителни променливи са последните стойности на всяка променлива.»

Автоматичната регресия работи при едно основно предположение. Че ценностите от близкото минало съдържат достатъчно информация за бъдещето. Но това може да не е вярно.

Можете да опитате да извлечете повече информация от последните данни. Например, текущата обобщена статистика е полезна при описване на скорошната динамика.

Автоматизирано проектиране на функции

Инженерингът на функции включва извличане и куриране на обяснителни променливи. Това е ключов етап във всеки проект за наука за данни. Качеството на характеристиките е централен аспект на представянето на моделите. Поради това учените по данни прекарват много време в този процес.

И все пак инженерството на функции често е ad hoc процес. Учените, занимаващи се с данни, създават функции въз основа на своите познания и опит в областта. Така че автоматизирането на част от този процес е желателно за практикуващите.

Нека видим как можете да направите това за многовариантни времеви редове.

Инженеринг на функции за многомерни времеви редове

Четене на данните

Ще използваме многовариантна времева серия, „събрана от интелигентен буй“ като казус [1]. Този буй е поставен на брега на Ирландия. Той улавя 9 променливи, свързани с океанските условия. Те включват „морска температура, височина на вълните и скорост на морската вода“, наред с други. Фигура 1 по-горе показва графика за първия месец на 2022 г.

Ето как можете да прочетете тези данни с помощта на pandas:

import pandas as pd

# skipping second row, setting time column as a datetime column
# dataset available here: https://github.com/vcerqueira/blog/tree/main/data
buoy = pd.read_csv('data/smart_buoy.csv', 
                   skiprows=[1], 
                   parse_dates=['time'])

# setting time as index
buoy.set_index('time', inplace=True)
# resampling to hourly data
buoy = buoy.resample('H').mean()
# simplifying column names
buoy.columns = [
    'PeakP', 'PeakD', 'Upcross',
    'SWH', 'SeaTemp', 'Hmax', 'THmax',
    'MCurDir', 'MCurSpd'
]

Фигура 1 по-горе показва графика за първия месец на 2022 г.

Целта е да се прогнозират бъдещите стойности на променливата SWH (значителна височина на вълната). Тази променлива често се използва за количествено определяне на височината на океанските вълни. „Един от случаите на използване на този проблем е да се оцени енергията, произведена от океанските вълни“. Този енергиен източник е все по-популярна алтернатива на невъзобновяемите енергийни източници.

Авторегресивен модел

Времевият ред е многомерен. Така че можете да използвате подход ARDL (автоматично регресивно разпределени закъснения), за да решите тази задача. Можете да научите повече за този метод в предишната ми публикация.

Ето как бихте приложили този метод.

import pandas as pd

from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_percentage_error as mape
from sklearn.multioutput import MultiOutputRegressor
from lightgbm import LGBMRegressor

# https://github.com/vcerqueira/blog/blob/main/src/tde.py
from src.tde import time_delay_embedding

target_var = 'SWH'

colnames = buoy.columns.tolist()

# create data set with lagged features using time delay embedding
buoy_ds = []
for col in buoy:
    col_df = time_delay_embedding(buoy[col], n_lags=24, horizon=12)
    buoy_ds.append(col_df)

# concatenating all variables
buoy_df = pd.concat(buoy_ds, axis=1).dropna()

# defining target (Y) and explanatory variables (X)
predictor_variables = buoy_df.columns.str.contains('\(t\-')
target_variables = buoy_df.columns.str.contains(f'{target_var}\(t\+')
X = buoy_df.iloc[:, predictor_variables]
Y = buoy_df.iloc[:, target_variables]

# train/test split
X_tr, X_ts, Y_tr, Y_ts = train_test_split(X, Y, test_size=0.3, shuffle=False)

# fitting a lgbm model without feature engineering
model_wo_fe = MultiOutputRegressor(LGBMRegressor())
model_wo_fe.fit(X_tr, Y_tr)

# getting forecasts for the test set
preds_wo_fe = model_wo_fe.predict(X_ts)

# computing the MAPE error
mape(Y_ts, preds_wo_fe)
# 0.238

Първо, времевият ред се трансформира в авторегресивен проблем. Това става с функцията time_delay_embedding. Прогностичната цел е да се прогнозират следващите 12 стойности на SWH (horizon=12). Обяснителните променливи са последните 24 стойности на всяка променлива в серията (n_lags=24).

LightGBM се обучава за всеки прогнозен хоризонт с помощта на директен подход. Директният метод е популярен подход за многоетапно прогнозиране напред. Той е внедрен в scikit-learn с името MultiOutputRegressor.

Кодът по-горе изгражда и тества авторегресивен модел. Обяснителните променливи включват само последните стойности на всяка променлива. Това води до средна абсолютна процентна грешка от 0,238. Нека да видим дали този резултат може да бъде подобрен с инженеринг на функции.

Това ръководство включва два подхода за извличане на характеристики от многовариантни времеви редове:

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

Извличане на едномерни характеристики

Можете да обобщите последните стойности на всяка променлива. Например изчисляване на пълзяща средна стойност, за да се обобщи скорошното ниво. Или подвижна дисперсия, за да добиете представа за скорошната степен на дисперсия.

import numpy as np

SUMMARY_STATS = {
    'mean': np.mean,
    'sdev': np.std,
}

univariate_features = {}
# for each column in the data
for col in colnames:
    # get lags for that column
    X_col = X.iloc[:, X.columns.str.startswith(col)]

    # for each summary stat
    for feat, func in SUMMARY_STATS.items():
        # compute that stat along the rows
        univariate_features[f'{col}_{feat}'] = X_col.apply(func, axis=1)

# concatenate features into a pd.DF
univariate_features_df = pd.concat(univariate_features, axis=1)

Може да искате да добавите още няколко статистики. Можете да направите това, като добавите функции към речника SUMMARY_STATS. Наличието на тези функции в един речник поддържа кода подреден.

Извличане на двумерни характеристики

Едномерната статистика пропуска потенциални взаимодействия между различни променливи. Можете да заснемете тази информация, като използвате процес на извличане на двумерни характеристики.

Идеята е да се изчислят функции за различни двойки променливи. И така, вие обобщавате общата динамика на тези двойки с двоична статистика.

Има два начина да направите това:

• Подвижни двоични статистики. Изчислете статистики, които приемат двойки променливи като вход. Например подвижната ковариация или подвижната корелация;
• Подвижна двоична трансформация, последвана от едномерна статистика. Трансформирайте двойка променливи в една променлива и обобщете тази променлива. Например изчисляване на кръстосаната корелация по елементи и след това вземане на нейната средна стойност.

Примери за подвижна двоична статистика включват ковариация, корелация или относителна ентропия.

Има много възможности за двоична трансформация. Например процентна разлика, кръстосана корелация или линейна конволюция между двойки променливи. След това тези трансформации се обобщават със статистически данни като средно или стандартно отклонение.

Ето скрипта, използван за извършване на тези два процеса на извличане наведнъж.

import itertools

import pandas as pd

from scipy.spatial.distance import jensenshannon
from scipy import signal
from scipy.special import rel_entr

from src.feature_extraction import covariance, co_integration

BIVARIATE_STATS = {
    'covariance': covariance,
    'co_integration': co_integration,
    'js_div': jensenshannon,
}

BIVARIATE_TRANSFORMATIONS = {
    'corr': signal.correlate,
    'conv': signal.convolve,
    'rel_entr': rel_entr,
}

# get all pairs of variables
col_combs = list(itertools.combinations(colnames, 2))

bivariate_features = []
# for each row
for i, _ in X.iterrows():
    # feature set in the i-th time-step
    feature_set_i = {}
    for col1, col2 in col_combs:
        # features for pair of columns col1, col2

        # getting the i-th instance for each column
        x1 = X.loc[i, X.columns.str.startswith(col1)]
        x2 = X.loc[i, X.columns.str.startswith(col2)]

        # compute each summary stat
        for feat, func in BIVARIATE_SUMMARY_STATS.items():
            feature_set_i[f'{col1}|{col2}_{feat}'] = func(x1, x2)

        # for each transformation
        for trans_f, t_func in BIVARIATE_TRANSFORMATIONS.items():

            # apply transformation
            xt = t_func(x1, x2)

            # compute summary stat
            for feat, s_func in SUMMARY_STATS.items():
                feature_set_i[f'{col1}|{col2}_{trans_f}_{feat}'] = s_func(xt)

    bivariate_features.append(feature_set_i)

bivariate_features_df = pd.DataFrame(bivariate_features, index=X.index)

Отново можете да добавите допълнителни трансформации или статистики. Това става чрез включването им в речниците BIVARIATE_TRANSFORMATIONS или BIVARIATE_STATS.

След като извлечете всички функции, вие ги свързвате към оригиналните обяснителни променливи. След това цикълът на обучение и тестване е като това, което сте правили преди.

# concatenating all features with lags
X_with_features = pd.concat([X, univariate_features_df, bivariate_features_df], axis=1)

# train/test split
X_tr, X_ts, Y_tr, Y_ts = train_test_split(X_with_features, Y, test_size=0.3, shuffle=False)

# fitting a lgbm model with feature engineering
model_w_fe = MultiOutputRegressor(LGBMRegressor())
model_w_fe.fit(X_tr, Y_tr)

# getting forecasts for the test set
preds_w_fe = model_w_fe.predict(X_ts)

# computing MAPE error
print(mape(Y_ts, preds_w_fe))
# 0.227

Това води до 0,227 средна абсолютна процентна грешка, което е подобрение. Подходът, който не използва инженеринг на функции, претърпя по-голяма загуба (0,238).

Избор на функция

Процесът на извличане по-горе води до общо 558 обяснителни променливи.

В зависимост от броя на променливите и обобщената статистика, това може да доведе до проблем с големи размери. Затова е важно да премахнете лошите или излишни функции от набора от данни.

Един от начините да направите това е да вземете най-важните функции и да повторите процеса на обучение с тях.

# getting the importance of each feature in each horizon
avg_imp = pd.DataFrame([x.feature_importances_
                        for x in model_w_fe.estimators_]).mean()

# getting the top 100 features
n_top_features = 100

importance_scores = pd.Series(dict(zip(X_tr.columns, avg_imp)))
top_features = importance_scores.sort_values(ascending=False)[:n_top_features]
top_features_nm = top_features.index

# subsetting training and testing sets by those features
X_tr_top = X_tr[top_features_nm]
X_ts_top = X_ts[top_features_nm]

# re-fitting the lgbm model
model_top_features = MultiOutputRegressor(LGBMRegressor())
model_top_features.fit(X_tr_top, Y_tr)

# getting forecasts for the test set
preds_top_feats = model_top_features.predict(X_ts_top)

# computing MAE error
mape(Y_ts, preds_top_feats)
# 0.229

Топ 100 функции водят до подобна производителност като пълните 558 функции.

Ето важността на първите 15 характеристики (други са пропуснати за краткост):

Най-важната характеристика е първият лаг на целевата променлива. Но някои от извлечените функции присъстват в този топ 15. Например, третата най-добра функция в SWH|Hmax_js_div. Това представлява отклонението на Jensen-Shannon между закъсненията на целевата променлива и закъсненията на Hmax. Петата най-добра характеристика е SeaTemp_sdev: стандартното отклонение на морската температура изостава. Ковариацията също е подходяща статистика за различни двойки променливи.

Друг начин за премахване на излишни функции е прилагането на корелационен филтър. Вие премахвате силно корелирани характеристики, за да намалите измерението на данните.

Обобщаване на целия времеви ред

Това ръководство се фокусира върху проблеми с прогнозиране с многовариантни времеви редове. И така, процесът на извличане на характеристики се прилага към много подпоследователности от времева серия. На всяка времева стъпка обобщавахте данните за последните 24 часа с набор от статистически данни.

Но можете също да приложите тези статистики, за да характеризирате цялата времева поредица наведнъж. Това може да е полезно, ако целта ви е да групирате набор от времеви редове. Първо, вие обобщавате всеки времеви ред с извличане на функции. След това прилагате алгоритъм за клъстериране към получените функции.

Ключови неща за вкъщи

• Многовариантното прогнозиране на времеви редове обикновено е авторегресивен процес
• Инженерингът на функции е ключова стъпка в проектите за наука за данни.
• Можете да подобрите многовариантни набори от данни от времеви серии с инженеринг на функции. Това включва изчисляване на едномерни и двумерни трансформации и обобщена статистика.
• Извличането на твърде много функции води до проблем с големи размери. Можете да използвате методи за избор на функции, за да премахнете нежелани функции.

Благодаря за четенето и ще се видим в следващата история!

Препратки

[1] Източник на данни: https://erddap.marine.ie/erddap/tabledap/IWaveBNetwork.html. Лиценз: Creative Commons Attribution 4.0 (https://erddap.marine.ie/erddap/info/IWaveBNetwork/index.html)

[2] Серкейра, Витор, Нуно Мониш и Карлос Соареш. „Жилетка: Автоматично проектиране на функции за прогнозиране.“ Машинно обучение (2021): 1–23.