Exibir código-fonte para Mínimos Quadrados

<!---
--->

Este o nome que se da ao ajuste ou ''fitting'' de uma função (polinômio) a um conjunto de dados.

Se <math>(X_i, Y_i)</math> com <math>i=1, N</math> representam o conjunto de dados (N) obtidos de um experimento (instrumento) ou 
de uma observação (por exemplo, em pesquisa de opinião ou censo) ou de uma simulação numérica.
E se suspeitamos que existe uma correlação entre os X (variável independente ou de entrada, controlada pelo experimento)
e os Y (cuja dependência com X queremos testar), primeiro colocamos os pontos num gráfico para ver se o conjunto forma
uma nuvem dispersa (quando não existe correlação aparente, isto é X e Y não conformam uma função), ou se existe 
correlação (os pontos parecem estar sobre alguma curva). <br>

=Equação linear=

[[Image:minquadlinear.png|right|frame|Exemplo de ajuste linear para um conjunto de pontos.]]

Sendo que um experimento foi realizado e temos <math>N</math> pontos, como descrito acima, e consideramos que um ajuste linear é coerente, uma reta deve ser construída para melhor representar estes pontos. Como mostrado na figura a baixo, para cada ponto, teremos um erro <math>\epsilon_i</math>, que é definido como a distância entre o ponto experimental e a curva (reta neste caso) teórica que desejamos ajustar, ou seja,

:<math>\epsilon_i = y_i - f(x_i)</math> ,

onde

<!---
:<math>f(x) = \alpha_0 + \alpha_1x + \alpha_2x^2 + \alpha_3x^3 + .... + \alpha_nx^n</math>
--->

:<math>f(x) = \alpha_0 + \alpha_1x  </math>

é a função que representa a curva de melhor ajuste. 

Para encontrar a reta que melhor se ajusta aos dados experimentais, desejamos minimizar o erro <math>\epsilon</math>. Como o erro pode ter tanto valores negativos quanto positivos, o que importa é minimizar o valor absoluto  de <math>\epsilon_i</math>. Isto poderia ser feito minimizando módulo de <math>\epsilon_i</math>, mas como a função módulo tem uma descontinuidade, é mais fácil minimizar o quadrado do erro. Para isto, definimos: 

:<math>S = \sum_{i=1}^N \epsilon_i^2</math>,
assim
:<math>S = \sum_{i=1}^N [y_i - f(x_i)]^2 = \sum_{i=1}^N [y_i - f(x_i;\alpha_0,\alpha_1)]^2 </math>.

Para obter a melhor reta que se ajusta aos dados experimentais, temos que minimizar <math>S</math> em relação às constantes da função <math>(\alpha_0,\alpha_1)</math>:

:<math>\frac{\partial S}{\partial \alpha_i} = 0</math>.

Como a reta possui apenas dois coeficientes, para o ajuste linear temos duas equações:

:<math>\frac{\partial S}{\partial \alpha_0} = \frac{\partial}{\partial \alpha_0}\sum_{i=1}^N [y_i - (\alpha_0 + \alpha_1x_i)]^2 = 0</math>
e
:<math>\frac{\partial S}{\partial \alpha_1} = \frac{\partial}{\partial \alpha_1}\sum_{i=1}^N [y_i - (\alpha_0 + \alpha_1x_i)]^2 = 0</math> .

Derivando as equações acima, temos que

:<math>\sum_{i=1}^N y_i - \sum_{i=1}^N \alpha_0 - \sum_{i=1}^N \alpha_1x_i = 0</math>
e
:<math>\sum_{i=1}^N y_ix_i - \sum_{i=1}^N \alpha_0x_i - \sum_{i=1}^N \alpha_1x_i^2 = 0</math> .

Assim,

:<math>\alpha_0\underbrace{\sum_{i=1}^N 1}_N + \alpha_1\underbrace{\sum_{i=1}^N x_i}_X = \underbrace{\sum_{i=1}^N y_i}_Y </math>
e
:<math>\alpha_0\underbrace{\sum_{i=1}^N x_i}_X + \alpha_1\underbrace{\sum_{i=1}^N x_i^2}_{X^2} = \underbrace{\sum_{i=1}^N y_ix_i}_{YX}</math> .

Lembre-se de que os valores <math>x_i</math> e <math>y_i</math> são conhecidos (são dados do problema). Desse modo, terminamos com um sistema linear para resolver, que na notação matricial fica

:<math>\begin{pmatrix}N & X \\ X & X^2\end{pmatrix} \begin{pmatrix} \alpha_0 \\ \alpha_1 \end{pmatrix} = \begin{pmatrix} Y \\ YX \end{pmatrix}</math> .


Cuidado com o fato que <math>(X^2 \neq X*X)</math> e <math>(YX \neq Y*X)</math>. Após construir a matriz, resolva com o método que mais lhe agrade (ha diversos métodos de solução de sistemas lineares, tais como a Regra de Cramer ou a eliminação Gaussiana).

=Equação quadrática=

[[Image:minquadquadratico.png|right|frame|Exemplo de ajuste quadrático para um conjunto de pontos.]]

Utilizando o mesmo método descrito para um ajuste linear, considerando que o melhor ajuste para um conjunto de pontos seja uma curva proveniente de função quadrática, temos que a função é dada por

:<math>f(x) = \alpha_0 + \alpha_1x + \alpha_2x^2</math> .

Desse modo, a soma do quadrado do erro fica

:<math>S = \sum_{i=1}^N \epsilon_i^2 = \sum_{i=1}^N [Y_i - (\alpha_0 + \alpha_1X_i + \alpha_2X_i^2)]^2</math>.

Após algumas contas, como feito na seção anterior, temos o sistema linear de 3 equações e 3 incógnitas para resolver:

:<math>\begin{pmatrix} N & X & X^2 \\ X & X^2 & X^3 \\ X^2 & X^3 & X^4 \end{pmatrix} \begin{pmatrix} \alpha_0 \\ \alpha_1 \\ \alpha_2 \end{pmatrix} = \begin{pmatrix} Y \\ YX \\ YX^2 \end{pmatrix} </math> .

Fique atento ao fato de que

:<math>X = \sum_{i=1}^N X_i\, , \, Y = \sum_{i=1}^N Y_i \, , \, X^2 = \sum_{i=1}^N X_i^2 \, , \, X^3 = \sum_{i=1}^N X_i^3 \, , \, X^4 = \sum_{i=1}^N X_i^4 \, , \, YX = \sum_{i=1}^N Y_iX_i \,\, e \,\, YX^2 = \sum_{i=1}^N Y_iX_i^2</math> .



=Polinômio de grau n=

Generalizando o procedimento acima, apresentado para polinômios de grau 1 e 2, podemos ajustar um conjunto de pontos com um polinômio de um grau específico <math>n</math>. Assim, a função será descrita por

:<math>f(x) = \alpha_0 + \alpha_1x + \alpha_2x^2 + \alpha_3x^3 + ... + \alpha_nx^n</math>

e a soma dos quadrados do erro é dada por

:<math>S = \sum_{i=1}^N \epsilon_i^2 = \sum_{i=1}^N [Y_i - f(X_i;\alpha_0,\alpha_1,...,\alpha_n)]^2</math>.

Ao final do procedimento, teremos um sistema linear de <math>n</math> equações e <math>n</math> incógnitas para resolver. O resultado deste sistema são os coeficientes :<math>\alpha_0, \alpha_1, \alpha_2.. \alpha_n</math> que compõem o polinômio que melhor se ajusta aos dados experimentais.

:<math>\begin{pmatrix}N & X & X^2 & \dots & X^n \\ X & X^2 & X^3 & \dots & X^{n+1} \\ X^2 & X^3 & X^4 & \dots & X^{n+2} \\ \vdots & \vdots & \vdots & \, & \vdots \\ X^n & X^{n+1} & X^{n+2} & \dots & X^{2n}\end{pmatrix} \begin{pmatrix} \alpha_0 \\ \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_n \end{pmatrix} = \begin{pmatrix} Y \\ YX \\ YX^2 \\ \vdots \\ YX^n \end{pmatrix}</math>



=Outros tipos de funções=

Dependendo do tipo de experimento, podem haver outras relações entre os pontos, como funções exponenciais.

==Exponencial 1==

Se os dados de um experimento se ajustarem bem  a uma função exponencial do tipo:

:<math>f(x) = \alpha_1 e^{-\alpha_2x} \, , \, \, \, \alpha_1\, , \, \alpha_2 > 0 </math> ,

definimos uma nova  função :

:<math>f_2(x) = \ln(f(x)) = \ln(\alpha_1 e^{-\alpha_2 x}) = \ln(\alpha_1) - \alpha_2 x</math>.

Assim, recaímos no problema do ajuste linear recém visto: 

:<math>f_2(x) = c_1 + c_2x</math>, com <math> c_1 = \ln(\alpha_1) </math> e  <math> c_2 = - \alpha_2</math>.

==Exponencial 2==

Se a função exponencial for do tipo:

:<math>f(x) = \alpha_1 \alpha_2^x</math>,

supondo <math>f(x)>0</math>, definimos:

:<math>f_2(x) = \ln(f(x)) = \ln(\alpha_1) + x\ln(\alpha_2)</math>.

Assim, como no caso anterior, voltamos para o problema de ajuste linear:

:<math>f_2(x) = c_1 + c_2x</math>,

com <math> c_1 = \ln(\alpha_1) </math> e  <math> c_2=\ln(\alpha_2)</math>.


==Algébrica==

Se a função for do tipo:

:<math>f(x) = \alpha_1 x^{\alpha_2}</math> ,

com <math>f(x)>0</math> e <math>x>0</math>, definimos:

:<math>f_2(x) = \ln(f(x)) = \ln(\alpha_1) + \alpha_2\ln(x)</math>.

e assim

:<math>f_2(x) = c_1 + c_2\ln(x)</math> ,

onde <math>c_1 = \ln(\alpha_1)</math> e <math>c_2 = \alpha_2</math>. Note também que os valores de x devem ser transformados em <math>\ln(x)</math> para ajustar os pontos.


----
=Código FORTRAN=
A seguir vemos uma possível implementação do método em linguagem F90.<br>
Observem a simplicidade do mesmo:

<pre>
 ! programa fortran para ajuste linear de conjunto de dados
 Implicit none
 Real :: xi,yi, x,y,xy,x2
 Real :: det,a,b

 n = 0;  x = 0;  y = 0;  xy = 0;  x2 = 0
 Do
    Read(*,*,end=100) xi,yi
    n = n + 1                          ! soma do numero de pontosd
    x  = x  + xi;      y =  y + yi     ! somatorio dos x e y
    x2 = x2 + xi**2;  xy = xy + xi*y   ! somatorio dos x**2 e x*y <- cuidado ha um erro aqui (compila mas ...
 End Do

 100 det = n*x2 - x**2
 a =  y*x2 - xy*x / det  ! <- outro erro aqui
 b = ...          / det  !    fica como exercicio

 print*, 'a=', a, 'b=', b
 end
</pre>
=Ajuste ponderado=
Dependendo da situação, convém fazer um ajuste levando em conta o erro associado a cada ponto, i.e., atribuindo maior peso para pontos com um erro baixo e menor peso para os pontos onde o erro é sabidamente maior. Ou seja, se definirmos <math>w_i</math> como o peso associado ao ponto <math>(X_i,Y_i)</math>, gostaríamos que ele seja maior quanto menor for o erro associado a este ponto.
Se <math>S_{y_i}</math> é o erro associado a este ponto, e considerando que o ajuste proposto é tal que minimiza a distância quadrática, podemos definir então <math>w_i</math> como:

<math>w_i=S_{y_i}^{-2}</math>

E o resíduo <math>\chi^2</math>, para o cálculo do ajuste ponderado, será dada por:

<math>\chi^2 = \sum_{i=1}^N (y_i - a - b x_i)^{2}w_i</math>

Aplicando o mesmo procedimento anterior para minimizar <math> \chi^2 </math>, obtemos as equações

:<math> \begin{bmatrix}
a[w] + b[xw]   & = &[yw] \\
a[xw] + b[x^2w] & = & [xyw] 
\end{bmatrix}</math>

E, portanto, os valores de <math>a</math> e <math>b</math> são:

<math>a = ([yw][x^2w] - [xyw][xw]) / \Delta</math>

<math>b = ([w][xyw]  -  [xw][yw]) / \Delta</math>

com <math>\Delta</math>:

<math>\Delta = [w][x^2w] - [xw]^2</math>

=Erro dos coeficientes=
Vimos como obter os coeficientes (a e b para uma reta) do ajuste de um conjunto de dados.<br /> 
Também como fazer esse ajuste quando os erros na variável dependente y não são todos iguais.<br />
Mas como saber se esses coeficientes são "bons". Ou seja, que margem de erro eles tem.<br />
Intuitivamente sabemos que quanto maior seja a dispersão dos <math>y_i</math> em volta
da curva do ajuste, maior será nossa incerteza sobre os coeficientes.

Vamos ver como traduzir isso de forma quantitativa. Voltando as expressões dos coeficientes a e b,
eles são funções de <math>x_i</math> e <math>y_i</math>, onde só os segundos são considerados como 
fonte de erro. Assim para ver como o erro neles propaga-se para os coeficientes, escrevemos:

<math>a = a(y_i) \Rightarrow \frac{\partial a}{\partial y_i} = \frac{1}{\Delta} \frac{\partial }{\partial y_i} \left\{[yw][x^2w] - [xyw][xw]\right\}</math>

<math>\frac{\partial a}{\partial y_i} = \frac{1}{\Delta}\left\{ w_i[x^2w] - x_i w_i [xw]\right\}</math>

Só os termos com y contribuem para a derivada. Como os <math>y_i</math> aparecem somados, ao derivar respeito do ''i''-esimo
sobra apenas o que multiplica ele

Para incluir o efeito do erro de cada <math>y_i</math> deveriamos somar i de 1 a N, mas como o erro pode ser para mais o menos
fazemos uma media quadrática deles:

<math>\Delta a = \sqrt{\sum_{i=1}^N \left( \frac{\partial a}{\partial y_i} \Delta y_i \right)^2}</math>

onde:
<math>\left(\frac{\partial a}{\partial y_i}\Delta y_i \right)^2 = \frac{1}{\Delta^2} \left\{ w_i^2[x^2w]^2 + x_i^2 w_i^2 [xw]^2 - 2w_i[x^2w] x_i w_i [xw]\right\}w_i^{-1}</math>

o somatório fica:

<math>\frac{1}{\Delta^2} \sum_{i=1}^N \left(w_i[x^2w]^2 + x_i^2 w_i [xw]^2 - 2 [x^2w] x_i w_i [xw]\right) = 
\frac{1}{\Delta^2} \left([w][x^2w]^2 + [x^2w][xw]^2 - 2 [x^2w][xw][xw]\right)</math>


e com mais algumas simplificações chegamos a simples relação:

<math>\Delta a = \sqrt{\frac{[x^2w]}{\Delta}}</math>

Analogamente para o b (que resulta ser mais fácil), se chega a:

<math>\Delta b = \sqrt{\frac{[w]}{\Delta}}</math>


Podemos interpretar essa expressões no caso sem ponderar, ou seja quando todos os erros são iguais:

<math>w=1/(\Delta y)^2</math>

:<math>\Rightarrow \Delta = w^2 (N[x^2] - [x]^2)= (w N \sigma)^2</math>

onde <math>\sigma^2 = <x^2> - <x>^2</math>

resultando:

<math>\Delta a = \frac{\Delta y \sqrt{<x^2>}}{\sigma \sqrt{N}}</math>

<math>\Delta b = \frac{\Delta y}{\sigma \sqrt{N}}</math>

Por tanto, três fatores determinam a qualidade do ajuste:

* E erro das medidas (<math>\Delta y</math>) que deve ser minimizado, porem está geralmente limitado pelo instrumento utililizado
* O número de medidas <math>N</math>, quanto maior, melhor, porem vemos que o erro dos coeficientes diminui com a raiz dele
* Por último, a dispersão da viariável dependente x (<math>\sigma</math>) também, quanto maior, melhor

Por último, consideremos o caso mais simples, ajuste sem ponderar e sem informação sobre <math>\Delta y</math>, apenas <math>N</math> pares de dados {<math>{x_i, y_i}</math>}.
Nesse caso podemos estimar o valor de  <math>\Delta y</math> pelo <math>\chi^2</math> assim:

 <math>\Delta y \approx \sqrt{\frac{\chi^2}{N-2}}</math>

Pois o resíduo <math>\chi^2</math> pode ser interpretado como o somatório dos erros de cada ponto, tomados como a distancia entre a medida e o ajuste.
O <math>N-2</math> em lugar de <math>N</math>, é pelo fato de que o ajuste já contem dois parâmetros (<math>a</math> e <math>b</math>
obtidos dos dados, então os erros individuais não são todos independentes.