##############################################################################
#                          FÍSICA COMPUTACIONAL II                           #
#                                   por                                      #
#                         Francisco  Carlos  Lavarda                         #
##############################################################################

Tutorial de MAXIMA: Ajuste de Curvas pelo Método dos Mínimos Quadrados
======================================================================

APARENTEMETE O MANUAL DE REFERÊNCIA EM PORTUGUÊS BRASILEIRO ESTÁ
DESATUALIZADO. DESTE MODO RECOMENDO O USO DA VERSÃO EM INGLÊS.


** 
Para os trabalhos com as funções para ajustes de curvas, carregue antes a
livraria lsquares:

load("lsquares");
 

** 
Defina a lista de pontos lista00 que pertencem à reta y=2x+1:

lista00:[[-1,-1],[-1/2,0],[0,1],[0.5,2],[1,3],[1.5,4]];


** 
Faça um gráfico a partir desta lista:

plot2d([discrete,lista00],[style,points]);

Para a versão wxMaxima também existe a função wxplot2d (que se usa da mesma
forma que plot2d).


** 
Faça a matriz matriz00 colocando os pontos de lista00 em colunas:

matriz00:matrix([-1,-1],[-1/2,0],[0,1],[0.5,2],[1,3],[1.5,4]);


** 
Calcule os coeficientes de reta que melhor se ajustem a estas colunas de
dados. Chamaremos a primeira coluna de x e a segunda de y. A curva a ser
ajustada é y=ax+b. Os coeficientes a serem calculados serão a e b.

lsquares_estimates(matriz00,[x,y],y=a*x+b,[a,b]);

Vemos que de os coeficientes calculados dão os valores esperados a=2 e
b=1.
  
Esta função lsquares_estimates permite que se ajuste uma curva qualquer a um
conjunto de pontos.

  
** 
Faça o gráfico dos pontos e da curva:

plot2d([[discrete,lista00],2*x+1],[x,-2,3],[y,-2,6],[style,points,lines]);


** 
Calcule o erro quadrático médio (eqm) para os coeficientes obtidos (2 e 1):

lsquares_residual_mse(matriz00,[x,y],y=a*x+b,[a=2,b=1]);


** 
Agora defina a lista lista01, que é uma alteração dos pontos de lista00,
de tal modo que tenhamos uma representação parecida com um experimento que
tem por lei teórica a equação y=2x+1:

lista01:[[-1.1,-0.9],[-0.4,0.1],[0.1,0.9],[0.4,2.1],[0.9,3.1],[1.6,3.9]];


** 
Faça o gráfico:

plot2d([discrete,lista01],[style,points]);


** 
Faça a matriz matriz01 colocando os pontos de lista01 em colunas:

matriz01:matrix([-1.1,-0.9],[-0.4,0.1],[0.1,0.9],[0.4,2.1],[0.9,3.1],[1.6,3.9]);


** 
Calcule os coeficientes de reta que melhor se ajustem a estas colunas de
dados. Chamaremos a primeira coluna de x e a segunda de y. A curva a ser
ajustada é y=ax+b. Os coeficientes a serem calculados serão a e b.

lsquares_estimates(matriz01,[x,y],y=a*x+b,[a,b]);

Vemos que os coeficientes calculados dão valores próximos a a=2 e b=1.
  

** 
Faça o gráfico dos pontos e da curva com os valores encontrados para a e
b:

plot2d([[discrete,lista01],(1722/907)*x+(28807/27210)],[x,-2,3],[y,-2,6],[style,points,lines]);


** 
Calcule o eqm para os coeficientes obtidos:

lsquares_residual_mse(matriz01,[x,y],y=a*x+b,[a=(1722/907),b=(28807/27210)]);

Veja que agora o eqm é diferente de zero.
  
   
** 
Não seria o caso, mas suponha que você suspeite que na realidade
uma curva parabólica traria uma descrição mais acurada para o fenômeno que
gerou este conjunto de pontos. Qual seria o eqm obtido nesta situação?
Veja você mesmo, com os comandos abaixo:

lsquares_estimates(matriz01,[x,y],y=a*x^2+b*x+c,[a,b,c]);
plot2d([[discrete,lista01],(-5/268)*x^2+(927527/486152)*x+(651217/607690)],[x,-2,3],[y,-2,6],[style,points,lines]);
lsquares_residual_mse(matriz01,[x,y],y=a*x^2+b*x+c,[a=(-5/268),b=(927527/486152),c=(651217/607690)]);

Veja que o eqm foi  menor, pois se permitiu que a curva de ajuste
pudessse se "flexionar" para minimizar o desvio.   


**
Agora vamos provocar um ajuste que nós sabemos de antemão ser muito ruim:
uma reta passando pela origem. Siga os comandos:

lsquares_estimates(matriz01,[x,y],y=a*x,[a]);
plot2d([[discrete,lista01],(1091/491)*x],[x,-2,3],[y,-2,6],[style,points,lines]);
lsquares_residual_mse(matriz01,[x,y],y=a*x,[a=(1091/491)]);

Observa-se um eqm muito maior que aquele obtido com uma equação de reta
mais geral.

Nota-se assim que o eqm serve como ferramenta para se decidir sobre a
qualidade do ajuste.


**
Até agora falamos de ajuste de curvas mais simples. Mas a função
lsquares_estimate permite o ajuste de qualquer tipo de curva. Por exemplo:

lsquares_estimates(matriz01,[x,y],(y+d)^2=a*x^2+b*x+c,[a,b,c,d]);
lsquares_estimates(matriz01,[x,y],y=a*exp(b*x),[a,b]);


**
Permite também o ajuste de curvas com mais de uma variável independente:

matriz02:matrix([1,1,1], [3/2,1,2], [9/4,2,1], [3,2,2], [2,2,1]);
lsquares_estimates(matriz02,[x,y,z],z=a*x+b*y+c,[a,b,c]);
lsquares_estimates(matriz02,[x,y,z],z=a*x+b*y+c*x*y+d,[a,b,c,d]);
lsquares_estimates(matriz02,[x,y,z],(z+e)^2=a*x+b*y+c*x*y+d,[a,b,c,d,e]);


**
Ao contrário da função lsquares_estimates que permite um ajuste para
qualquer tipo de curva, a função plsquares busca obrigatoriamente o melhor
ajuste para uma função do tipo polinômio. Ou seja, você não tem idéia do
tipo de curva que melhor descreve o conjunto de pontos que você tem em mãos
e você quer somente o melhor polinômio possível para aquele conjunto.

Para os trabalhos a seguir, carregue a livraria plsquares:

load("plsquares");

Comecemos por perguntar: qual teria sido o melhor polinômio de primeiro grau
(ou seja, uma reta) para o conjunto de pontos descritos por
lista00/matriz00?

plsquares(matriz00,[x,y],y,1);

A resposta é exatamente o resultado que esperávamos. O "determination
coefficient" varia entre 0 e 1: quanto mais próximo de 1, melhor é o ajuste
do polinômio encontrado.

Como um teste para a função plsquares, façamos a pergunta: se eu permitir
que o polinômio a ser obtido seja de grau 2, o que se obteria?

plsquares(matriz00,[x,y],y,2);

Percebe-se que o resultado é o mesmo obtido anteriormente: o melhor
polinômio é uma reta mesmo, em vista do fato que os pontos dados estão
localizados exatamente em cima de uma reta.

Dependendo da versão do software, o uso de plsquares para este conjunto de
pontos pedindo por um polinômio de grau>1 pode apresentar um erro; muitas
vezes ele pode ser solucionado com a seguinte sequência de comandos:
(1) debugmode(true); (para entrar no modo de depuração) 
(2) plsquares(matriz00,[x,y],y,2); (depois de acusar o erro entra em modo 
    de debug)
(3) :continue; (para continuar o cálculo e apresentar o resultado) 
(4):quit; (para sair do modo de depuração).

Agora: qual o melhor polinômio de grau=1 para lista01/matriz01?

plsquares(matriz01,[x,y],y,1);

Faça os cálculos e veja que o resultado é levemente diferente do que o
obtido com lsquares_estimates em que se força um polinômio ax+b.

Já para os dados em matriz02 os resultados são idênticos aos obtidos com
lsquares_estimates:

plsquares(matriz02,[x,y,z],z,1);

IMPORTANTE: qual é o grau máximo do polinômio pedido? O grau máximo está
limitado ao número de linhas da matriz, menos 1 (n-1).

**
Vejamos um exemplo mais elaborado para plsquares com a matriz03:

matriz03:matrix([-2,2],[-1,2.1],[0,2.5],[1,5.7],[2,29.3]);

que, para efeito de aprendizado, posso dizer que são pontos próximos à curva
y=2+(1/2)exp(2x).

Podemos procurar por um ajuste exatamente deste tipo, que é y=a+b.exp(cx)
através da função lsquares_estimate:

lsquares_estimates(matriz03,[x,y],y=a+b*exp(c*x),[a,b,c]);

que dará como resposta valores de a, b e c muito próximos àqueles da curva
da qual tomamos os pontos em sua proximidade.  Agora calcula-se o mse:

lsquares_residual_mse(matriz03,[x,y],y=a+b*exp(c*x),[a=2.00835193603933,b=0.49908109633254,c=2.000783158067087])

Agora tentemos o melhor polinômio possível para este conjunto de pontos:

plsquares(matriz03,[x,y],y,4);

e calculemos o mse desta curva com:

lsquares_residual_mse(matriz03,[x,y],y=(295*x^4+804*x^3+401*x^2+60*(x+19.76000000000001))/480,[a=1,b=1]);

Vemos que ambos os erros são da ordem de 10**(-4); mesmo assim, analisando
os gráficos para as duas funções e os pontos, você acha que é possível
escolher qual é o melhor ajuste?

**
Nota sobre o uso da função lsquares_residual_mse:
Quando temos uma função já conhecida e para a qual queremos calcular o eqm
em relação a um conjunto de pontos, aparentemente temos um problema, pois
não existem valores de parâmetros a serem repassados para a função. Assim,
para calcular o eqm da função y=3x+2 em relação aos pontos em "matriz",
podemos adotar um dos dois procedimentos abaixo:

lsquares_residual_mse(matriz,[x,y],y=a*x+b,[a=3,b=2]);

em que re-criamos uma função desconhecida y=ax+b e damos os parâmetros que já
eram conhecidos, ou

lsquares_residual_mse(matriz,[x,y],y=3*x+2,[a=1,b=1]);

em que os valores para a e b podem ser quaisquer que não interferem com o
cálculo.