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 $(X_{i},Y_{i})$ com $i=1,N$ 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).

Equação linear

Exemplo de ajuste linear para um conjunto de pontos.

Sendo que um experimento foi realizado e temos $N$ 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 $\epsilon _{i}$, que é definido como a distância entre o ponto experimental e a curva (reta neste caso) teórica que desejamos ajustar, ou seja,

$\epsilon _{i}=Y_{i}-f(X_{i})$ ,

onde

$f(x)=\alpha _{0}+\alpha _{1}x$

é 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 $\epsilon$. Como o erro pode ter tanto valores negativos quanto positivos, o que importa é minimizar o valor absoluto de $\epsilon _{i}$. Isto poderia ser feito minimizando módulo de $\epsilon _{i}$, mas como a função módulo tem uma descontinuidade, é mais fácil minimizar o quadrado do erro. Para isto, definimos:

$S=\sum _{i=1}^{N}\epsilon _{i}^{2}$,

assim

$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}$.

Para obter a melhor reta que se ajusta aos dados experimentais, temos que minimizar $S$ em relação às constantes da função $(\alpha _{0},\alpha _{1})$:

${\frac {\partial S}{\partial \alpha _{i}}}=0$.

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

${\frac {\partial S}{\partial \alpha _{0}}}={\frac {\partial }{\partial \alpha _{0}}}\sum _{i=1}^{N}[Y_{i}-(\alpha _{0}+\alpha _{1}X_{i})]^{2}=0$

e

${\frac {\partial S}{\partial \alpha _{1}}}={\frac {\partial }{\partial \alpha _{1}}}\sum _{i=1}^{N}[Y_{i}-(\alpha _{0}+\alpha _{1}X_{i})]^{2}=0$ .

Derivando as equações acima, temos que

$\sum _{i=1}^{N}Y_{i}-\sum _{i=1}^{N}\alpha _{0}-\sum _{i=1}^{N}\alpha _{1}X_{i}=0$

e

$\sum _{i=1}^{N}Y_{i}X_{i}-\sum _{i=1}^{N}\alpha _{0}X_{i}-\sum _{i=1}^{N}\alpha _{1}X_{i}^{2}=0$ .

Assim,

$\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}$

e

$\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_{i}X_{i}} _{YX}$ .

Lembre-se de que os valores $X_{i}$ e $Y_{i}$ são conhecidos (são dados do problema). Desse modo, terminamos com um sistema linear para resolver, que na notação matricial fica

${\begin{pmatrix}N&X\\X&X^{2}\end{pmatrix}}{\begin{pmatrix}\alpha _{0}\\\alpha _{1}\end{pmatrix}}={\begin{pmatrix}Y\\YX\end{pmatrix}}$ .

Cuidado com o fato que $(X^{2}\neq X*X)$ e $(YX\neq Y*X)$. 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

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

$f(x)=\alpha _{0}+\alpha _{1}x+\alpha _{2}x^{2}$ .

Desse modo, a soma do quadrado do erro fica

$S=\sum _{i=1}^{N}\epsilon _{i}^{2}=\sum _{i=1}^{N}[Y_{i}-(\alpha _{0}+\alpha _{1}X_{i}+\alpha _{2}X_{i}^{2})]^{2}$.

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

${\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}}$ .

Fique atento ao fato de que

$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_{i}X_{i}\,\,e\,\,YX^{2}=\sum _{i=1}^{N}Y_{i}X_{i}^{2}$ .

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 $n$. Assim, a função será descrita por

$f(x)=\alpha _{0}+\alpha _{1}x+\alpha _{2}x^{2}+\alpha _{3}x^{3}+...+\alpha _{n}x^{n}$

e a soma dos quadrados do erro é dada por

$S=\sum _{i=1}^{N}\epsilon _{i}^{2}=\sum _{i=1}^{N}[Y_{i}-f(X_{i};\alpha _{0},\alpha _{1},...,\alpha _{n})]^{2}$.

Ao final do procedimento, teremos um sistema linear de $n$ equações e $n$ incógnitas para resolver. O resultado deste sistema são os coeficientes :$\alpha _{0},\alpha _{1},\alpha _{2}..\alpha _{n}$ que compõem o polinômio que melhor se ajusta aos dados experimentais.

${\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}}$

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:

$f(x)=\alpha _{1}e^{-\alpha _{2}x}\,,\,\,\,\alpha _{1}\,,\,\alpha _{2}>0$ ,

definimos uma nova função :

$f_{2}(x)=\ln(f(x))=\ln(\alpha _{1}e^{-\alpha _{2}x})=\ln(\alpha _{1})-\alpha _{2}x$.

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

$f_{2}(x)=c_{1}+c_{2}x$, com $c_{1}=\ln(\alpha _{1})$ e $c_{2}=-\alpha _{2}$.

Exponencial 2

Se a função exponencial for do tipo:

$f(x)=\alpha _{1}\alpha _{2}^{x}$,

supondo $f(x)>0$, definimos:

$f_{2}(x)=\ln(f(x))=\ln(\alpha _{1})+x\ln(\alpha _{2})$.

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

$f_{2}(x)=c_{1}+c_{2}x$,

com $c_{1}=\ln(\alpha _{1})$ e $c_{2}=\ln(\alpha _{2})$.

Algébrica

Se a função for do tipo:

$f(x)=\alpha _{1}x^{\alpha _{2}}$ ,

com $f(x)>0$ e $x>0$, definimos:

$f_{2}(x)=\ln(f(x))=\ln(\alpha _{1})+\alpha _{2}\ln(x)$.

e assim

$f_{2}(x)=c_{1}+c_{2}\ln(x)$ ,

onde $c_{1}=\ln(\alpha _{1})$ e $c_{2}=\alpha _{2}$. Note também que os valores de x devem ser transformados em $\ln(x)$ para ajustar os pontos.

---

Código FORTRAN

A seguir vemos uma possível implementação do método em linguagem F90.
Observem a simplicidade do mesmo:

 ! 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

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 $w_{i}$ como o peso associado ao ponto $(X_{i},Y_{i})$, gostaríamos que ele seja maior quanto menor for o erro associado a este ponto. Se $S_{y_{i}}$ é o erro associado a este ponto, e considerando que o ajuste proposto é tal que minimiza a distância quadrática, podemos definir então $w_{i}$ como:

$w_{i}=S_{y_{i}}^{-2}$

E a variável $\chi$, para o cálculo do ajuste ponderado, será dada por:

$\chi =\sum _{i=1}^{N}(Y_{i}-a-bX_{i})^{2}w_{i}$

Aplicando o mesmo procedimento anterior para minimizar $\chi$, obtemos as equações

${\begin{bmatrix}a[w]+b[Xw]&=&[Yw]\\a[Xw]+b[X^{2}w]&=&[XYw]\end{bmatrix}}$

E, portanto, os valores de $a$ e $b$ são:

$\Delta =[w][X^{2}w]-[Xw]^{2}$

$a=[Yw][X^{2}w]-[XYw][Xw]/\Delta$

$b=[w][XYw]-[Xw][Yw]/\Delta$

Erro dos coeficientes

Vimos como obter os coeficientes (a e b para uma reta) do ajuste de um conjunto de dados.
Também como fazer esse ajuste quando os erros na variável dependente y não são todos iguais.
Mas como saber se esses coeficientes são "bons". Ou seja, que margem de erro eles tem.
Intuitivamente sabemos que quanto maior seja a dispersão dos $y_{i}$ 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 $x_{i}$ e $y_{i}$, onde só os segundos são considerados como fonte de erro. Assim para ver como o erro neles propaga-se para os coeficientes, escrevemos:

$a=a(y_{i})\Rightarrow {\frac {\partial a}{\partial y_{i}}}={\frac {1}{\Delta }}{\frac {\partial }{\partial y_{i}}}\left\{[Yw][X^{2}w]-[XYw][Xw]\right\}$

${\frac {\partial a}{\partial y_{i}}}={\frac {1}{\Delta }}\left\{w_{i}[X^{2}w]-x_{i}w_{i}[Xw]\right\}$

pois só os termos com y contribuem para a derivada. e como os $y_{i}$ aparecem somados, ao derivar respeito do i-esimo sobra apenas o que multiplica ele

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

$\Delta a={\sqrt {\sum _{i=1}^{N}\left({\frac {\partial a}{\partial y_{i}}}\Delta y_{i}\right)^{2}}}$

onde: $\left({\frac {\partial a}{\partial y_{i}}}\Delta y_{i}\right)^{2}={\frac {1}{\Delta ^{2}}}\left\{w_{i}^{2}[X^{2}w]^{2}+x_{i}^{2}w_{i}^{2}[Xw]^{2}-2w_{i}[X^{2}w]x_{i}w_{i}[Xw]\right\}w_{i}^{-1}$

o somatório fica:

${\frac {1}{\Delta ^{2}}}\sum _{i=1}^{N}\left(w_{i}[X^{2}w]^{2}+x_{i}^{2}w_{i}[Xw]^{2}-2[X^{2}w]x_{i}w_{i}[Xw]\right)={\frac {1}{\Delta ^{2}}}\left([w][X^{2}w]^{2}+[X^{2}w][Xw]^{2}-2[X^{2}w][Xw][Xw]\right)$

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

$\Delta a={\sqrt {\frac {[X^{2}w]}{\Delta }}}$

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

$\Delta b={\sqrt {\frac {[w]}{\Delta }}}$

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

$w=1/(\Delta y)^{2}$

$\Rightarrow \Delta =w^{2}(N[X^{2}]-[X]^{2})=(wN\sigma )^{2}$

onde $\sigma ^{2}=<x^{2}>-<x>^{2}$

resultando:

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

$\Delta b=\Delta y{\sqrt {\frac {1}{N\sigma ^{2}}}}$

Ou seja quanto maior o número e mais espalhadas no eixo X as medidas melhor