Um breve tutorial sobre Séries Temporais Nebulosas — Parte II (com estudo de caso em Energia Solar)

Na primeira parte desse tutorial eu falei o que eram séries temporais, conjuntos nebulosos e por fim o que eram as Séries Temporais Nebulosas — STN, com uma breve introdução da biblioteca pyFTS. Como em todo estudo introdutório, utilizamos os métodos mais simples com dados simples, com o objetivo de facilitar o entendimento da idéia geral da abordagem.

Mas agora, partindo do pressuposto que os leitores já possuem o background necessário, podemos avançar um pouco e brincar com algo mais útil. Primeiro vamos rever intuitivamente o conceito de viés e variância, na sequência vamos nos aprofundar nas STN, entendendo os seus principais hyperparâmetros e seu impacto na acurácia dos modelos e por fim vamos fazer um estudo de caso em séries temporais de radiação solar, usadas para prever a geração de energia fotovoltailca.

Todos os exemplos de código desse tutorial estão disponíveis no Google Colab, no endereço http://bit.ly/tutorial_fts_colab2 . Sinta-se a vontade para entrar em contato, tirar dúvidas e dar sugestões! E agora vamos começar, né?

O sinal, o ruído, o viés e a variância

O treinamento de modelos de inteligência computacional é um embate para separar o sinal do ruído, o comportamento geral das especificidades locais e erros aleatórios. Como em todo estimador baseado em inteligência computacional, há um embate entre o viés e a variância, entre o sub-ajuste e o sobre ajuste.

Vamos supor uma série temporal Y, e seus valores individuais y(t)∈Y, e um estimador de Y, uma função ŷ(t+1) = f( y(t) ). O objetivo de f é predizer Y da melhor maneira possível fazendo com que a diferença ε(t) entre o valor real y(t) e o valor predito ŷ(t) tenda a zero, em outras palavras ε(t) = y(t) — ŷ(t).

Intuitivamente, um estimador enviesado/viciado (biased) é aquele cujo valor esperado de ε(t) é diferente de zero — E[ε(t)] ≠ 0. O viés é um “desvio” sistemático do valor correto. Lembre que o estimador frequentemente terá desvios do valor real, e isso é esperado, mas na média esses desvios devem tender a zero. O enviesamento é típico do sub ajuste do modelo (underfitting), quando o modelo não foi capaz de aprender o sinal, o componente da série temporal que nos interessa.

Fonte

A variância por outro lado tem a ver com o poder de generalização do modelo, em especial com dados que não lhe foram apresentados durante o treinamento. A alta variância está relacionada com o sobre ajuste do modelo (overfitting), quando o modelo começa a aprender especificidades da amostra de treino que não generalizam bem para os dados de teste. Em suma: o modelo aprendeu o ruído dos dados.

Sabe-se que não é possível ter um modelo não-enviesado e de baixavariância, e o ajuste ideal (best fit) é chegar a um equilíbrio entre esses termos — e esse é o desafio dos modelos de estimação.

Parâmetros das STN

Fonte

Diversos parâmetros determinam o melhor ajuste de um modelo de STN, mas principalmente o particionamento e a ordem. Esses dois parâmetros sozinhos respondem por 90% (valor empírico) da acurácia do modelo.

1. Particionamento

O particionamento é composto por três parâmetros, aqui listados segundo sua importância:

1a) Número de partições (ou conjuntos nebulosos)

Simplesmente o parâmetro mais influente na acurácia do modelo. Quanto mais conjuntos nebulosos mais precisa é a captação das características da série temporal. E a armadilha reside justamente aqui:

  • Poucos conjuntos geram sub ajuste, simplificando demasiadamente o sinal;
  • Conjuntos em excesso geram sobre ajustes, capturando até o ruído dos dados;
Diversos particionamentos para a função seno
Acurácia dos diversos particionamentos para a função seno

A quantidade de conjuntos é um parâmetro que deve ser testado. Em séries que foram diferenciadas 10 é um bom número para começar. Nos demais caso 35 é um bom número para começar.

1b) Tipo de particionamento

Existem muitos tipos de particionamento, desde o particionamento em Grade (GridPartitioner) onde todos os conjuntos tem o mesmo formato, passando pelos particionadores onde conjuntos tem tamanhos distintos — como os baseados em entropia e clusterização. Não vou me aprofundar nessa discussão aqui, mas para os curiosos há vários exemplos dos tipos de particionamento no repositório PYFTS/notebooks.

Sempre comece com o particionamento em Grid, depois se for o caso, explore os demais tipos.

1c) Função de pertinência

Esse é um parâmetro que tem pouca influência real na acurácia do modelo, mas dependendo do caso você pode ter uma justificativa para usar uma função gaussiana ou trapezoidal no lugar da função triangular, que é a padrão.

Uma das justificativas pode ser a quantidade de parâmetros (a gaussiana usa 2, a triangular 3 e a trapezoidal 4), a legibilidade do modelo ou mesmo outras questões ligadas à natureza do processo e dos dados.

Novamente não vou me aprofundar nessa discussão aqui, dê uma olhada no repositório PYFTS/notebooks.

2. Ordem

A ordem do modelo é o segundo parâmetro mais importante das STN, dado que elas são modelos autoregressivos (usam valores passados para prever os próximos). A ordem é o tamanho da memória do modelo, quanta informação passada é necessária para descrever acontecimentos futuros.

Para escolher esses parâmetros é importante estar familiarizado com o conceito de Função de Autocorrelação — ACF. A ACF não apenas é capaz de indicar a ordem como também os índices dos lags mais importantes.

2a) Quantidade de lags

A ordem é a quantidade de lags (valores passados) que são usados pelos modelos. Isso é realmente muito importante e depende do tipo de série temporal que está sendo modelada. A pergunta aqui é: quantos lags eu preciso para que o modelo aprenda os padrões temporais, ciclos, sazonalidades, etc. Olhe o ACF e veja quantos são os lags mais significativos.

Acurácia dos modelos segundo sua ordem (nº de lags)

Porém, quanto mais lags maior o modelo (principalmente se o número de partições for grande!), e quanto maior o modelo mais lento será o aprendizado e inferência.

Pela minha experiência, dificilmente é necessário mais do que 3 lags para se descrever uma série temporal. Mas é claro: tudo depende dos dados.

2b) Escolha dos índices do lags

Por padrão os últimos lags, em ordem, são escolhidos. Mas dependendo das sazonalidades da série temporal esses podem não ser os melhores lags, então olhe o ACF e veja quais os índices lags mais significativos.

3. Tipos de métodos

A literatura acerca dos métodos de STN é muito variada, mas duas diferenciações são extremamente importantes:

3a) Com Pesos x Sem Pesos

Os pesos aumentam a acurária do modelo, ponderando quais os conjuntos nas regras do modelo são mais influentes. Eles diminuem um pouco a legibilidade do modelo, mas nada demasiado. Se tiver que escolher, sempre prefira os modelos com pesos!

No exemplo abaixo podemos comparar o modelo HOFTS (sem pesos), WHOFTS (com pesos no consequente da regras) e o PWFTS (com pesos no consequente e no precedente das regras):

3b) Monovariados x Multivariados

Na maior parte do tempo só dispomos de séries temporais com uma variável — a variável endógena. Outras vezes essa variável vem auxialada de outras informações (as variáveis exógenas) das quais podemos tirar proveito.

A data, geralmente associada às medições da série temporal, é uma variável preciosíssima no caso de dados sazonais — como dados sociais, ambientais, etc.

A primeira coisa é saber se há de fato correlação entre as variáveis, então utilize uma matriz de correlação para verificar. As correlções apontam relações lineares simples então não deve ser as únicas ferramentas que você deva utilizar, a entropia cruzada é uma boa alternativa.

Outra sugestão: Se você tem uma série monovariada você pode enriquecer seu modelo criando uma série multivariada a partir dela, onde as outras variáveis são transformações da variável exógena. Por exemplo, você pode ter uma série multivariada com a variável original (endógena) e a variável endógena diferenciada, provendo uma informação extra sobre as flutuações recentes dos valores.

Estudo de Caso: Radiação Solar

Vamos pôr a mão na massa agora né? Só para relembrar: todos os códigos e resultados estão disponíveis no endereço http://bit.ly/tutorial_fts_colab2 . Vou usar os dados de radiação solar da base SONDA — Sistema de Organização Nacional de Dados Ambientais (INPE), especificamente da estação de Brasília/DF. Os dados de radiação estão na variável ‘glo_avg’, amostrados por minuto, 24 horas por dia, de 2013 à 2014 para o conjunto de treino e o ano de 2015 para teste.

Como é esperado, os dados são um pouco ruidosos e não há grande necessidade da amostragem por minuto. Fiz uma limpesa nos dados e retirei as outras variáveis (cuja correlação com nossa variável principal não é tão relevante) e deixei apenas a data e a variável endógena. A etapa seguinte foi reduzir o volume de dados e a limpeza dos ruídos, eu fiz isso amostrando a média horária da série. Os dados já preprocessados estão no endereço https://data.world/petroniocandido/sonda-bsb-averaged-hourly .

Intuitivamente a radiação solar tem dois ciclos principais, o ciclo diário (movimento de rotação) e o ciclo anual (o movimento de translação). A maior incerteza que afeta a predição da radiação são as nuvens!

Amostras do ciclo anual
Amostras do ciclo diário
Função de autocorrelação para os primeiros 24 lags, correspondento às 24 horas anteriores

Métodos para todos os gostos e bolsos

Como vocês já sabem os métodos na biblioteca pyFTS ficam no pacote models. Vamos analisar então como diversos modelos de STN se saem na modelagem dessa série temporal. No mundo real agora faríamos uma otimização dos hiperparâmetros do modelo testando inúmeras combinações de todos os parâmetros citados acima. Na pyFTS, usando o pacote hyperparam, você pode optar por uma GridSearch distribuída, que fará um produto cartesiano de todas os valores possíveis de cada hiperparâmetro e testar cada um deles. Você também poderia utilizar um algoritmo genético (hyperparam.Evolutive) ou a biblioteca hyperopt. Mas vamos simplificar, ok?

Vamos particionar nossa variável endógena em grade (GridPartitioner) usando 35 partições, separados em 5 subgrupos, MB — Muito Baixo, B — Baixo, M — Médio, A — Alto e MA — Muito Alto, cada subgrupo com 7 níveis. Essa nomenclatura tornará o modelo mais interpretável gerando regras como:

B5 , B6 → B6 , M0 , M1

E esta pode ser lida como:

SE y(t-1) é Baixo (subnível 5) E y(t) é Baixo (subnível 6) ENTÃO y(t+1) será Baixo(subnível 6) OU Médio (subnível 0) OU Médio (subnível 1).

Particionamento em grade da série temporal de radiação solar

Métodos Monovariados

Já temos o particionamento então vamos explorar os métodos. Optamos pelos métodos de alta ordem (ordem > 1) com e sem pesos e todos eles foram testados com a ordem de 1 à 3. Os métodos e exemplos das regras geradas seguem abaixo:

  • hofts.HighOrderFTS: Método de alta ordem, sem pesos
B4,MB0 → MB0,MB1
B5,MB0 → MB0,MB1
  • hofts.WeightedHighOrderFTS: Método de alta ordem, com pesos no consequente das regras
B4,MB0 → MB0 (0.5), MB1 (0.5)
B5,MB0 → MB0 (0.5), MB1 (0.5)
  • pwfts.ProbabilisticWeightedFTS: Método de alta ordem, com pesos probabilísticos no precedente e consequente das regras
0.003) B0,MB1 → (0.876)MB0, (0.103)MB1, (0.015)MB2, (0.006)MB3, (0.001)MB4
(0.003) B0,MB2 → (0.003)B0, (0.003)B1, (0.003)B2, (0.0)B3, (0.0)B4, (0.787)MB0, (0.164)MB1, (0.03)MB2, (0.002)MB3, (0.002)MB4, (0.005)MB5, (0.001)MB6
Amostra do desempenho dos modelos univariados

Métodos Multivariados

Na pyFTS os modelos multivariados usam Pandas Dataframes como entrada de dados (os modelos monovariados usam lists ou arrays Numpy) e os métodos estão no pacote models.multivariate. Cada variável tem seu próprio particionamento e outros parâmetos. Deve-se primeiro criar um objeto do tipo Variable, informando o nome da variável, a coluna de dados no dataframe e o particionador. Em seguida adicionam-se as variáveis exógenas e a endógena no modelo usando a função append_variable. A variável endógena deve ser setada na propriedade target_variable.

from pyFTS.models.multivariate import common, variable, mvfts
from pyFTS.models.seasonal import partitioner as seasonal
from pyFTS.models.seasonal.common import DateTime
sp = {'seasonality': DateTime.day_of_year , 'names': ['Jan','Fev','Mar','Abr','Mai','Jun','Jul', 'Ago','Set','Out','Nov','Dez']}
mes = variable.Variable("Month", data_label="data", partitioner=seasonal.TimeGridPartitioner, npart=12,data=train_mv, partitioner_specific=sp)
radiacao = variable.Variable("Radiation", data_label="glo_avg", alias='rad',partitioner=Grid.GridPartitioner, npart=35, data=train_mv)
model = mvfts.MultivariateFTS()
model.append_variable(mes)
model.append_variable(radiacao)
model.target_variable = radiacao
model.fit(train_mv)

Vamos utilizar só o campo data e dele extrair duas variáveis exógenas: o mês e a hora. O particionamento das variáveis segue abaixo:

Você deve estar se perguntando: Porquê os conjuntos de horas e meses estão sobrepostos? Afinal essa é uma informação muito exata! Resposta: Isso é lógica nebulosa minha gente! Partimos do princípio que coisas próximas se influenciam e têm comportamentos parecidos, daí as regras de uma hora específica também sofrem influência das horas vizinhas e meses idem.

Vamos agora dar uma olhada nos métodos e ver exemplos das regras geradas por eles:

  • mvfts.MultivariateFTS: Método de primeira ordem, sem pesos;
Jan,8hs,MB0 → MB2,MB3,MB4,MB1,MB5
Jan,8hs,MB1 → MB2,MB3,B1,MB6,MB5,MB4,B0,MB1,B2,B4
  • wmvfts.WeightedMultivariatedFTS: Método de primeira ordem, com pesos;
Jan,8hs,MB0  → MB2 (0.353), MB1 (0.253), MB4 (0.147), MB3 (0.207), MB5 (0.04)
Jan,8hs,MB1  → MB2 (0.276), MB3 (0.172), MB1 (0.198), MB5 (0.083), MB4 (0.151), MB6 (0.021), B0 (0.036), B4 (0.005), B1 (0.036), B2 (0.021)
  • cmvfts.ClusteredFTS: Método de alta ordem, com pesos;
Jan11hsM3,Jan12hsB1  → Jan13hsMB6 (1.0)
Jan12hsB1,Jan13hsMB6 → Jan15hsMB3 (1.0)

Comparando os modelos

Obviamente o principal critério para avaliar um modelo preditivo é sua acurácia. Vamos utilizar a Raiz do Erro Médio Quadrático (Root Mean Squared Error — RMSE) e quanto menor esse valor melhor. Mas outro critério importante é a parcimônia: os modelos mais simples serão sempre mais desejáveis que os modelos complexos. Então temos como objetivo minimizar o RMSE e o tamanho do modelo. Vamos aos cinco modelos com melhor desempenho:

Era de se esperar que os modelos multivariados tivessem um desempenho melhor, afinal dispõem de mais informação. Mas veja que o modelo mais parcimonioso é o PWFTS!

As possibilidades não param aqui. Poderíamos agora resolver fazer um ensemble dos melhores modelos (models.ensemble), ou usar conjuntos nebulosos não estacionários (models.nonstationary). Mas isso fica para outra oportunidade!

Ampliando o horizonte de predição

Quando utilizamos modelos preditivos também estamos interessados em não apenas prever o próximo valor imediato mas uma sequência desses valores. Isso é bem legal pois, com apenas uns poucos lags de informação podemos gerar uma estimativa para muitos períodos à frente! Para o caso específico da previsão de radiação solar vamos fixar o nosso horizonte de predição em 48 horas.

Para um modelo monovariado isso é simples, basta utilizar o parâmetro steps_ahead da função predict. Esse parâmetro retroalimentará os valores de saída na entrada, pelas próximas steps_ahead iterações.

predicoes = model.predict(dados, steps_ahead=48)

Para modelos multivariados há uma complicação a mais: não estamos gerando somente valores da variável endógena mas das exógenas também. Como isso depende da natureza de cada variável a função predict exigirá, além do parâmetro steps_ahead, o parâmetro generators, que é um dicionário que deverá conter uma chave para cada variável exógena cujos valores são funções que recebem o valor atual da variável e devem retornar o próximo valor. No nosso exemplo isso é bem simples dado que os valores das variáveis exógenas são datetimes e basta adicionar uma hora ao valor de entrada:

generator = lambda x : x + pd.to_timedelta(1, unit='h')
predicoes = model.predict(dados, steps_ahead=48, 
generators={'Hour': generator, 'Month': generator})
Previsões para as próximas 48 horas, usando os dois melhores modelos

Concluindo

Agora você já tem conhecimento suficiente em STN e na biblioteca pyFTS para aplicá-los nos seus projetos pessoais. É claro que ainda falta falar de muitas coisas: previsão por intervalo e probabilística, transformações de dados, auditoria dos modelos, modelos incrementais e não estacionários… Mas isso fica para os próximas posts!

Não hesite em entrar em contato se precisar de ajuda, tirar dúvidas, tiver sugestões ou críticas. Até a próxima pessoal!