A previsão precisa de eventos futuros é uma tarefa promissora e interessante em muitas áreas: da previsão do tempo à fintech (cotações de ações, taxas de câmbio). Hoje, o aprendizado de máquina pode reduzir significativamente o tempo e o trabalho envolvidos na tomada de decisões de gerenciamento. Por cerca de meio ano,nossa equipe de Ciência de Dados da NORBIT experimentou vários modelos de aprendizado de máquina para resolver problemas de classificação e regressão e otimizar processos de negócios na esfera b2b. Mas quando a tarefa de prever a série temporal apareceu, os materiais disponíveis sobre esse tópico na rede não foram suficientes para desenvolver uma solução rápida.A essência da tarefa era prever com máxima precisão (erro absoluto médio em porcentagem, ou MAPE <3%) o consumo horário de uma cidade inteira pelos próximos três dias (72 horas) para as necessidades de uma grande empresa geradora de eletricidade.O erro máximo permitido de 3% é um indicador muito alto para a precisão das previsões. Antes do uso do aprendizado de máquina, o erro estava na faixa de 3-23%, o que levava diretamente a perdas financeiras, uma vez que, com a geração insuficiente de energia elétrica, era necessário adquirir capacidades adicionais no mercado atacadista de eletricidade (que é muito mais caro que a auto-gerado), com superprodução - vendido no mercado mais barato do que eles poderiam vender para a cidade. Além disso, o efeito negativo do desvio da carga planejada da carga real é uma alteração na frequência da corrente alternada na rede.Muitos fatores impedem a obtenção de alta precisão de previsão, por exemplo, com um aumento repentino de temperatura, as pessoas imediatamente procuram consoles de aparelhos de ar condicionado e, quando caem, recebem aquecedores nos mezaninos; durante as férias ou grandes jogos de futebol incluem TVs, etc. Enquanto alguns eventos são cíclicos e potencialmente previsíveis, outros estão completamente ausentes. Suponhamos que, de repente, devido a um acidente, a siderúrgica pare de funcionar e todas as previsões se transformem em abóbora (em nossa cidade-alvo não havia grandes empresas de produção). Havia um exemplo bastante interessante de um problema ao fazer previsões no setor de energia, quando um dos picos nas mudanças na quantidade de eletricidade gerada foi causado por um navio encalhado em um lago, e o capitão concordou com os proprietários das usinas hidrelétricas para reduzir a descarga de água.Naturalmente, o modelo de aprendizado de máquina não pôde prever esse evento com base em dados históricos.A boa notícia para nós é que o cliente forneceu esses "dias especiais" nos termos de referência que permitiram um erro médio de previsão de 8%.Tudo o que o cliente forneceu são os dados históricos do consumo horário de eletricidade por 3 anos e o nome da cidade.A primeira tarefa para preparar a construção de um modelo de aprendizado de máquina é a visualização e a busca de fatores adicionais que podem influenciar a previsão. O grau de tal influência é convenientemente avaliado visualmente usando um mapa de calor da correlação de atributos. O quadrado escuro, neste caso, significa a relação inversa incondicional das quantidades (quanto maior um valor, menor o outro e vice-versa), branco - direto:import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
sns.heatmap(df.corr());
Tanto quanto eu sei, nos modelos usados para prever o consumo de eletricidade na Federação Russa, dois fatores são usados - a história do consumo e a previsão de temperatura. No inverno e no outono, o consumo de energia é mais alto - temperatura baixa e horas de luz do dia relativamente curtas. O fator “dia da semana” também afeta a quantidade de energia gasta - nos fins de semana, diminui acentuadamente. Os picos de consumo ocorrem na quarta ou quinta-feira.No período de um dia da semana, os valores máximos de consumo de energia ocorrem de 11 a 12 horas e caem gradualmente com uma queda acentuada após as nove da noite. 
Tudo é como nos livros didáticos:Programações diárias de consumo diário de inverno e verão para centros industriais e culturais. FonteModelagem
Profeta
A primeira idéia de como tornar a tarefa o mais rápida possível é tirar o Profeta do Facebook. Eu já tinha experiência em usá-lo e lembrei-o como "algo que funciona de forma simples e rápida". O Profeta deseja ver as colunas "ds" e "y" no quadro de dados dos pandas, ou seja, data no formato AAAA-MM-DD HH: MM: SS e a meta é consumo por hora, respectivamente (exemplos de trabalho estão na documentação ).df = pd.read_pickle('full.pkl')
pjme = df.loc[:'2018-06-01 23:00'].rename(columns={'hour_value': 'y'})
pjme['ds'] = pjme.index
split_date = '2018-03-01 23:00'
pjme_train = pjme.loc[pjme.index <= split_date]
pjme_test = pjme.loc[pjme.index > split_date]
playoffs = pd.DataFrame({
'holiday': 'playoff',
'ds': df[(df.rest_day == 1)|(df.is_weekend == 1)].index,
'lower_window': 0,
'upper_window': 1,
})
model = Prophet(holidays=playoffs)
model.fit(pjme_train.reset_index())
forecast = model.predict(df=pjme_test.reset_index())
As previsões do "profeta" pareciam objetivas, mas o erro era 7 a 17% maior do que o permitido, portanto, não realizei mais experimentos nessa estrutura.Sarimaax
A segunda tentativa de resolver esse problema foi usar o SARIMAX. O método consome bastante tempo para selecionar coeficientes, mas foi possível reduzir o erro de previsão em comparação com o Profeta para 6-11%. Infelizmente, apenas o resultado das previsões semanais permaneceu, mas foram esses dados preditivos que usei como recursos adicionais nos modelos de impulso.Para a previsão, você precisa selecionar os parâmetros SARIMA (p, d, q) (P, D, Q, s):- p é a ordem autoregressiva (AR), que permite adicionar os valores passados da série temporal (pode ser expresso como "amanhã, provavelmente nevará se nevou nos últimos dias");
- d é a ordem de integração dos dados de origem (determina o número de pontos de tempo anteriores a serem subtraídos do valor atual);
- q é a ordem da média móvel (MA), que permite definir o erro do modelo na forma de uma combinação linear dos valores de erro observados anteriormente;
- P - p, mas sazonal;
- D - d, mas sazonal;
- Q - q, mas sazonal;
- s é a dimensão da sazonalidade (mês, trimestre, etc.).
Na expansão sazonal do histórico semanal, você pode ver a tendência e a sazonalidade, o que permite que o seasonal_decompose seja exibido:import statsmodels.api as sm
sm.tsa.seasonal_decompose(df_week).plot()
De acordo com os gráficos de autocorrelação e autocorrelação parcial, ele selecionou as aproximações iniciais dos parâmetros para o modelo SARIMAX: p = 0, P = 1, q = 2, Q = 2.sm.graphics.tsa.plot_acf(df_week.diff_week[52:].values.squeeze(), lags=60, ax=ax)
sm.graphics.tsa.plot_pacf(df_week.diff_week[52:].values.squeeze(), lags=60, ax=ax)
E os valores finais selecionados são (0, 2, 2) x (1, 2, 0, 52). Como resultado, o gráfico de previsão de consumo ficou assim:Impulsionar
Nesse ponto, ficou claro que, sem o uso de fatores externos adicionais, a precisão necessária não poderia ser alcançada. Primeiro de tudo, esses são fatores climáticos: temperatura, pressão, umidade, força e direção do vento, precipitação. A primeira tentativa foi uma tentativa de comprar um arquivo do tempo em Yandex. Os caras têm excelente API e suporte técnico, mas especificamente em nossa cidade havia muitas lacunas nos dados. Como resultado, comprei o arquivo no serviço openweathermap.org por US $ 10. É importante notar que, apesar do clima ser um fator muito útil, os erros de previsão obviamente estragam o modelo final do MAPE. Encontrei informações de que a solução correta para esse problema estaria em treinamento para usar não os valores históricos reais do tempo, mas previsões do tempo por três dias, mas, infelizmente, não tive essa oportunidade. Também adicionei sinais da hora do dia, dia da semana, feriados, número de série do dia do ano, além de um grande número de atrasos e valores médios para períodos anteriores.Todos os dados após o dimensionamento (MinMaxScale) e One Hot Encoding (os valores de cada item categórico se tornam colunas separadas com 0 e 1, portanto, o item com três valores se torna três colunas diferentes e a unidade estará em apenas um deles) que usei em um pequeno competição de três modelos populares de impulso XGBoost, LightGBM e CatBoost. XGBoost
Muito de bom foi escrito sobre o XGBoost. Apresentado na conferência SIGKDD em 2016, ele fez um splash e ainda mantém uma posição de liderança. Material interessante para quem nunca ouviu falar.Vamos preparar os dados para o modelo: from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn.model_selection import train_test_split
def gen_data(
stop_day = '2018-06-01 23:00',
drop_columns = [],
test_size=0.15,
is_dummies=True,
is_target_scale=False):
df = pd.read_pickle('full.pkl')
df = df.loc[:stop_day]
scaler = MinMaxScaler()
target_scaler = StandardScaler()
y = df.hour_total
if is_target_scale:
y = target_scaler.fit_transform(y.values.reshape(-1, 1)).reshape(1,-1)[0]
X = df.drop(['hour_value', *drop_columns], axis=1)
num_columns = X.select_dtypes(include=['float64']).columns
X[num_columns] = scaler.fit_transform(X[num_columns])
if is_dummies:
X = pd.get_dummies(X)
train_count_hours = len(X) - 72
valid_count_hours = len(X) - int(len(X) * 0.2)
X_test = X[train_count_hours:]
y_test = y[train_count_hours:]
X = X[:train_count_hours]
y = y[:train_count_hours]
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=test_size, random_state=42)
train_dates = X_train.index
val_dates = X_test.index
return X, y, X_train, X_val, y_train, y_val, X_test, y_test, target_scaler
Estou ensinando um modelo usando os melhores hiperparâmetros selecionados anteriormente:X, y, X_train, X_val, y_train, y_val, X_test, y_test, _ = gen_data(
stop_day = stop_day,
drop_columns = drop_columns,
test_size=0.1,
is_dump=True,
is_target_scale=False)
params_last = {
'base_score': 0.5,
'booster': 'gbtree',
'colsample_bylevel': 1,
'colsample_bynode': 1,
'colsample_bytree': 0.4,
'gamma': 0,
'max_depth': 2,
'min_child_weight': 5,
'reg_alpha': 0,
'reg_lambda': 1,
'seed': 38,
'subsample': 0.7,
'verbosity': 1,
'learning_rate':0.01
}
reg = xgb.XGBRegressor(**params_last, n_estimators=2000)
print(X_train.columns)
reg.fit(X_train, y_train,
eval_set=[(X_val, y_val)],
early_stopping_rounds=20,
verbose=False)
y_pred = reg.predict(X_test)
Lightgbm
O LightGBM é uma estrutura da Microsoft cuja principal vantagem é a velocidade de aprendizado em conjuntos de dados realmente grandes. E também, ao contrário do XGBoost, o LightGBM pode trabalhar com categorias, usa menos memória. Aqui então há mais.import lightgbm as lgb
from lightgbm import LGBMRegressor
from sklearn.metrics import mean_squared_error
stop_day = '2018-06-01 23:00'
start_day_for_iterate = datetime.strptime(stop_day, '%Y-%m-%d %H:%M')
test_size = 0.2
X, y, X_train, X_val, y_train, y_val, X_test, y_test, _ = gen_data(
stop_day = stop_day,
drop_columns = drop_columns,
test_size=test_size,
is_dump=True,
drop_weekend=False,
is_target_scale=False)
model = LGBMRegressor(boosting_type='gbdt',
num_leaves=83,
max_depth=-1,
learning_rate=0.008,
n_estimators=15000,
max_bin=255,
subsample_for_bin=50000,
min_split_gain=0,
min_child_weight=3,
min_child_samples=10,
subsample=0.3,
subsample_freq=1,
colsample_bytree=0.5,
reg_alpha=0.1,
reg_lambda=0,
seed=38,
silent=False,
nthread=-1)
history = model.fit(X_train, y_train,
eval_metric='rmse',
eval_set=[(X_val, y_val)],
early_stopping_rounds=40,
verbose = 0)
y_pred = model.predict(X_test, num_iteration=model.best_iteration_)
Catboost
O CatBoost é uma biblioteca avançada de aumento de gradiente em árvores de decisão de código aberto "da Yandex . A diferença vantajosa do modelo é a conveniência de trabalhar com conjuntos de dados que contêm atributos categóricos - você pode transferir dados que contêm categorias para o modelo sem conversão para números e revela padrões por conta própria nada precisa ser distorcido, e a qualidade da previsão permanece alta.cat_features=np.where(X.dtypes == 'category')[0]
eval_dataset = Pool(X_test, y_test)
model = CatBoostRegressor(learning_rate=0.5,
eval_metric='RMSE',
leaf_estimation_iterations=3,
depth=3)
model.fit(X_train, y_train,
eval_set=(X_val, y_val),
cat_features=cat_features,
plot=True,
verbose=True)
y_pred = model.predict(X_test)
XGBoost vs. LightGBM vs. Catboost
Para não repetir vários artigos, darei uma tabela comparativa a partir daqui .Para calcular o erro MAPE, tirei o último mês conhecido (28 dias) e, movendo a janela com uma resolução de uma hora, fiz uma previsão para as próximas 72 horas.E na minha competição improvisada os prêmios foram:3º lugar: meu favorito - XGBoost - com a menor qualidade de previsão;2º lugar: CatBoost como o mais conveniente e "pronto para uso";1º lugar: LightGBM como o mais rápido e preciso.Para uma comparação mais precisa dos modelos, usei R2 (R-quadrado ou coeficiente de determinação, mostrando o quanto a variação condicional do modelo difere da variação de valores reais) e RMSLE (Root Logarithmic Square Squared Mean Root, ou root mean square logarithmic Error, ou o erro logarítmico quadrado médio da raiz, que é, de fato, a distância entre dois pontos no plano - valor real e previsto). Mas tudo isso, é claro, em relação à minha tarefa. Em outros dados em mãos diferentes, tudo pode ser completamente diferente. Da próxima vez, vou compartilhar outra história sobre prever o fluxo de funcionários para a tarefa de planejar o recrutamento em nossa empresa.A propósito, nossa equipe está crescendo e estamos procurando pessoas!