Um exemplo real
Passo a passo para o ajuste de uma curva real
Vamos lá?
Já vimos um exemplo trabalhado e construímos a teoria. Nesta aula vamos aplicar tudo isto a um exemplo real.
Neste exemplo, vamos buscar uma curva que descreva bem o comportamento do nível dos reservatórios do sistema Cantareira, um dos principais reservatórios de água que abastecem a região metropolitana de São Paulo.
Os dados foram obtidos a partir do Portal dos Mananciais da Sabesp, ao longo de vários anos. Neste portal, as informações são atualizadas diariamente. Em minha página publico um acompanhamento desses dados, atualizado regularmente. Abaixo, reproduzo a imagem gerada no dia 01/12/2020.
![](img/cantareira.png)
Usando a informação desse portal, foi criado um arquivo de texto contendo leituras desde 2004 até o dezembro de 2020. Cada leitura está em uma linha e contém quatro colunas numéricas:
- Dia da leitura (valor inteiro entre 1 e 31)
- Mês da leitura (valor inteiro entre 1 e 12)
- Ano da leitura (valor inteiro entre 2004 e 2020)
- Percentual ocupado do volume do reservatório (valor real positivo)
Faça o download do arquivo de dados para poder acompanhar o desenvolvimento e fazer seus próprios experimentos.
Carregando e entendendo os dados
A primeira tarefa é inspecionar o conjunto de dados e conseguir exibi-lo. Para isso, vamos iniciar carregando os dados com o Octave.
load "cantareira.dat" whos Variables in the current scope: Attr Name Size Bytes Class ==== ==== ==== ===== ===== C 2708x4 86656 double Total is 10832 elements using 86656 bytes
Com isso já descobrimos que o dado é composto de 2708 leituras. Logo exibi-las em uma tabela é impraticável e pouco informativo. Seria muito mais útil criar um gráfico.
plot(C(:,4))
Esse gráfico foi construído sem explicitar qual vetor representa as ordenadas. Neste caso, o Octave utiliza com ordenada o índice da amostra, 1 a 2708. Já podemos tirar algumas conclusões comparando este gráfico com o exibido no topo da página. Claramente, as primeira 500 amostras devem ser anteriores a 2015. As amostras seguintes (aproximadamente 2200) correspondem a um período menor (5 anos). Portanto, a amostragem do dado não é regular, ou seja as amostras contidas no dado não foram coletadas a intervalos regulares de tempo. No início da série histórica temos poucas amostras, mais espaçadas no tempo, e a partir de 2015 a amostragem deve ficou mais densa. Para exibir o gráfico corretamente precisamos utilizar como vetor de ordenadas a informação correta da data de registro de cada amostra.
Para ajustarmos uma função aos dados precisamos representar a variável dependente (porcentagem do volume do reservatório, o vetor
C(1:4,:)
ans =
1.0000 9.0000 2004.0000 51.6504
15.0000 9.0000 2004.0000 49.5617
1.0000 10.0000 2004.0000 47.2409
15.0000 10.0000 2004.0000 46.0805
Vemos que a primeira leitura foi feita em 01/09/2004. Esse seria o dia 1. A leitura seguinte foi no feita em 15/09/2004, o que corresponderia ao dia 15. Já a terceira e quarta leituras corresponderiam aos dias corridos 31 e 46, respectivamente. Observe também que
C(end,:)
ans =
1.0000 12.0000 2020.0000 47.3183
A última leitura foi tomada em 01/12/2020, que em termos de dias corridos corresponde ao dia 5935. Converter datas em dias corridos a partir de uma data inicial é algo chato de fazer, por conta da diferença de dias em meses e anos bissextos. Além disso, trabalhar com esse intervalo grande para
Ao invés de pensar em dias corridos, podemos indexar as leituras em termos do ano contínuo, ou seja, uma amostra tomada em 01/01/2010 corresponderia a
x = ( C(:,3) + ((C(:,2)-1)*(365.25/12)+C(:,1))/365.25 ) - 2012; [min(x) max(x)] ans = -7.3306 8.9194 y = C(:,4); plot(x+2012,y)
Agora o vetor
Apesar de totalmente negligenciada nos textos clássicos sobre Cálculo Numérico, quando trabalhamos com dados reais, esta tarefa de carregar, inspecionar e preparar os dados é fundamental para a obtenção de um bom resultado.
Resolução do problema
Com os dados preparados, passamos ao problema de ajuste, em si. O primeiro a fazer é a escolha do espaço de aproximação, ou seja, das funções
A função
Usando como base para o espaço de aproximação
m = 2708; P = [ones(m,1) sin(2*pi*x) cos(2*pi*x)]; c = P\y c = 50.3609 5.3851 -3.9387 phi = P*c; plot(x+2012,y,'k', x+2012,phi,'r')
Claramente, conseguimos captar tanto o valor médio (por volta de 50), quanto as oscilações anuais, mas ainda estamos longe de ajustar bem a função. Certamente precisaremos de mais graus de liberdade para ajustar o comportamento complexo observado. É difícil antever quais exatamente serão as funções de base que vão dar origem a um bom ajuste. Neste ponto vamos nos permitir alguma liberdade em acrescentar mais termos, considerando bases
N = 5; P = ones(m,2*N+1); for k=1:N P(:,[2*k 2*k+1]) = [sin(2*k*pi*x) cos(2*k*pi*x)]; endfor c = P\y; phi = P*c; plot(x,y,x,phi)
Ao aumentar o quantidade de termos do tipo
whos P Variables in the current scope: Attr Name Size Bytes Class ==== ==== ==== ===== ===== P 2708x11 238304 double Total is 29788 elements using 238304 bytes P(:,12:13) = [sin(x) cos(x)]; c = P\y; phi = P*c; plot(x,y,x,phi)
Melhorou. Vamos adicionar alguns monômios, na esperança de que isto ajude a representar esse comportamento global. Neste último teste consideramos
Neste exemplo, fizemos considerações sobre período de oscilação observado no dado e usamos isso para guiar a escolha das funções de base. Talvez fosse mais adequado utilizar outra ferramenta como a transformada de Fourier discreta, que naturalmente lida com a decomposição do dado em funções trigonométricas de diferentes períodos. Mas isto não seria imediato neste problema, por conta da amostragem não regular. Além disso, uma vantagem da técnica de ajuste é que podemos mesclar funções distintas. O resultado da última figura foi obtido combinando-se 16 funções, entre monômios e funções trigonométricas. Será que apenas 16 termos da transformada de Fourier discreta produziriam resultado semelhante?
Em resumo, o problema de ajuste de curvas reside em dois aspectos fundamentais, compreender e preparar seu dado, e fazer boas escolhas para a base do espaço de aproximação.
1. Faça o download do dado usado nessa aula e tente reproduzir o que foi feito. Depois, experimente outras bases para o espaço de aproximação. Tente ajustar uma curva que represente apenas o comportamento médio dos dados, sem se preocupar em descrever oscilações anuais.
2. Tente determinar uma função que ajuste bem o comportamento da cotação do dólar. Faça o download do arquivo de dados. Neste arquivo há o registro da cotação do dólar desde 01/02/1999, quando entrou em vigor o sistema de câmbio flutuante no Brasil, até meados de outubro de 2021, momento em que o país amargava uma situação lamentável.