Mínimos Quadrados

Em 1809, Carl Friedrich Gauss (1777-1855) publicou um artigo no Werke, 4, 1-93, demonstrando que a melhor maneira de determinar um parâmetro desconhecido de uma equação de condições é minimizando a soma dos quadrados dos resíduos, mais tarde chamado de Mínimos Quadrados por Adrien-Marie Legendre (1752-1833). Em abril de 1810, Pierre-Simon Laplace (1749-1827) apresenta no memoir da Academia de Paris, ["Mémoire sur les approximations des formules qui sont fonctions de très-grands nombres, et sur leur application aux probabilités (suite)". Mémoires l'Institut 1809 (1810), 353-415, 559-565. Oeuvres 12 p.301-345, p.349-353] a generalização a problemas com vários parâmetros desconhecidos.

Um programa de mínimos quadrados sempre começa com a minimização da soma:

$ {S \equiv \sum_{i=1}^N (y^o_i - y_i)^2}$ (1)
onde chamamos de
yoi = valores observados de y
yi = valores calculados de y
ou seja, mínimos quadrados implica em minimizar os quadrados dos resíduos. Estes resíduos podem ser normalizados pelo número de pontos ou pela soma dos pesos empregados.
Residuos
Por que este critério é considerado um bom critério e não simplismente minimizar os resíduos ou o cubo dos resíduos? A resposta formal é que os mínimos quadrados são corretos se os resíduos tiverem uma distribuição gaussiana (normal).
Gaussiana
Distribuição gaussiana (normal) para uma variável x, com média m e desvio padrão s.
É simples notar que se minimizarmos os resíduos diretamente, um grande resíduo negativo pode ser anulado por um grande resíduo positivo, enquanto que com o quadrado minimizamos os módulos das diferenças.

Suponhamos que temos um conjunto de dados y com uma distribuição normal:

P(y) = $ {\frac{{1}}{{\sigma \sqrt{2\pi}}}}$exp$ [{-\frac{1}{2\sigma^2} (y-\bar{y})^2}.$ - $ {\frac{{1}}{{2\sigma^2}}}$$ ({y-\bar{y}}.$y - $ \bar{{y}}$$ {y-\bar{y}})^{2}_{}$$ {-\frac{1}{2\sigma^2} (y-\bar{y})^2}]$
onde
P(y) = probabilidade de obter o valor y  
y = quantidade a ser observada  
$ \bar{{y}}$ = valor médio de y  
$ \sigma$ = desvio padrão de y  

Por exemplo,
a probabilidade de se encontrar uma medida entre -$ \sigma$ e +$ \sigma$ é de 68,3%,
a probabilidade de se encontrar uma medida entre -1,5$ \sigma$ e +1,5$ \sigma$ é de 86,6%,
a probabilidade de se encontrar uma medida entre -2$ \sigma$ e +2$ \sigma$ é de 95,4%,
a probabilidade de se encontrar uma medida entre -2,5$ \sigma$ e +2,5$ \sigma$ é de 98,76%,
a probabilidade de se encontrar uma medida entre -3$ \sigma$ e +3$ \sigma$ é de 99,74%.

Suponhamos que medimos o valor de y várias vezes, obtendo uma série de valores {yi}. A probabilidade de observar este conjunto é dada por

P$ ({y_1,y_2,y_3,\ldots,y_N}.$y1, y2, y3,..., yN$ .{y_1,y_2,y_3,\ldots,y_N})$ = $ \{{\frac{1}{\sigma \sqrt{2\pi}}\exp [-\frac{1}{2\sigma^2}
(y_1-\bar{y})^2]}.$$ {\frac{{1}}{{\sigma \sqrt{2\pi}}}}$exp$ [{-\frac{1}{2\sigma^2}
(y_1-\bar{y})^2}.$ - $ {\frac{{1}}{{2\sigma^2}}}$$ ({y_1-\bar{y}}.$y1 - $ \bar{{y}}$$ .{y_1-\bar{y}})^{2}_{}$$ .{-\frac{1}{2\sigma^2}
(y_1-\bar{y})^2}]$$ .{\frac{1}{\sigma \sqrt{2\pi}}\exp [-\frac{1}{2\sigma^2}
(y_1-\bar{y})^2]}\}$$ \{{\frac{1}{\sigma \sqrt{2\pi}}\exp [-\frac{1}{2\sigma^2}
(y_2-\bar{y})^2]}.$$ {\frac{{1}}{{\sigma \sqrt{2\pi}}}}$exp$ [{-\frac{1}{2\sigma^2}
(y_2-\bar{y})^2}.$ - $ {\frac{{1}}{{2\sigma^2}}}$$ ({y_2-\bar{y}}.$y2 - $ \bar{{y}}$$ .{y_2-\bar{y}})^{2}_{}$$ {-\frac{1}{2\sigma^2}
(y_2-\bar{y})^2}]$$ {\frac{1}{\sigma \sqrt{2\pi}}\exp [-\frac{1}{2\sigma^2}
(y_2-\bar{y})^2]}\}$...  
    ...$ \{{\frac{1}{\sigma \sqrt{2\pi}}\exp [-\frac{1}{2\sigma^2}
(y_N-\bar{y})^2]}.$$ {\frac{{1}}{{\sigma \sqrt{2\pi}}}}$exp$ [{-\frac{1}{2\sigma^2}
(y_N-\bar{y})^2}$ - $ {\frac{{1}}{{2\sigma^2}}}$$ ({y_N-\bar{y}}.$yN - $ \bar{{y}}$$ .{y_N-\bar{y}})^{2}_{}$$ .{-\frac{1}{2\sigma^2}
(y_N-\bar{y})^2}]$$ {\frac{1}{\sigma \sqrt{2\pi}}\exp [-\frac{1}{2\sigma^2}
(y_N-\bar{y})^2]}\}$  
  = $ ({\frac{1}{\sigma \sqrt{2\pi}}}.$$ {\frac{{1}}{{\sigma \sqrt{2\pi}}}}$$ .{\frac{1}{\sigma \sqrt{2\pi}}})^{N}_{}$exp$ [{-\frac{1}{2\sigma^2}
\sum_{i=1}^N (y_i-\bar{y})^2}.$ - $ {\frac{{1}}{{2\sigma^2}}}$$ \sum_{{i=1}}^{N}$$ ({y_i-\bar{y}}.$yi - $ \bar{{y}}$$ .{y_i-\bar{y}})^{2}_{}$$ .{-\frac{1}{2\sigma^2}
\sum_{i=1}^N (y_i-\bar{y})^2}]$  
Queremos agora saber qual é o melhor valor de $ \bar{{y}}$, o valor médio. O melhor valor será aquele que maximiza a probabilidade de obter os valores observados {yi}, ou seja, o melhor valor de $ \bar{{y}}$ é obtido colocando a derivada da probabilidade como nula:
$ {\frac{{d}}{{d\bar{y}}}}$$ [{P(y_1,y_2,y_3,\ldots,y_N)}.$P$ ({y_1,y_2,y_3,\ldots,y_N}.$y1, y2, y3,..., yN$ .{y_1,y_2,y_3,\ldots,y_N})$$ .{P(y_1,y_2,y_3,\ldots,y_N)}]$ = 0

Ou seja

$ {\frac{{d}}{{d\bar{y}}}}$$ \{{(\frac{1}{\sigma \sqrt{2\pi}})^N \exp [-\frac{1}{2\sigma^2} \sum_{i=1}^N (y_i-\bar{y})^2]}.$$ ({\frac{1}{\sigma \sqrt{2\pi}}}.$$ {\frac{{1}}{{\sigma \sqrt{2\pi}}}}$$ .{\frac{1}{\sigma \sqrt{2\pi}}})^{N}_{}$exp$ [{-\frac{1}{2\sigma^2}
\sum_{i=1}^N (y_i-\bar{y})^2}.$ - $ {\frac{{1}}{{2\sigma^2}}}$$ \sum_{{i=1}}^{N}$$ ({y_i-\bar{y}}.$yi - $ \bar{{y}}$$ .{y_i-\bar{y}})^{2}_{}$$ .{-\frac{1}{2\sigma^2}
\sum_{i=1}^N (y_i-\bar{y})^2}]$$ .{(\frac{1}{\sigma \sqrt{2\pi}})^N \exp [-\frac{1}{2\sigma^2} \sum_{i=1}^N (y_i-\bar{y})^2]}\}$ = 0

$ ({\frac{1}{\sigma \sqrt{2\pi}}}.$$ {\frac{{1}}{{\sigma \sqrt{2\pi}}}}$$ .{\frac{1}{\sigma \sqrt{2\pi}}})^{N}_{}$exp$ [{-\frac{1}{2\sigma^2}
\sum_{i=1}^N (y_i-\bar{y})^2}.$ - $ {\frac{{1}}{{2\sigma^2}}}$$ \sum_{{i=1}}^{N}$$ ({y_i-\bar{y}}.$yi - $ \bar{{y}}$$ .{y_i-\bar{y}})^{2}_{}$$ .{-\frac{1}{2\sigma^2}
\sum_{i=1}^N (y_i-\bar{y})^2}]$$ {\frac{{d}}{{d\bar{y}}}}$$ [{-\frac{1}{2\sigma^2}
\sum_{i=1}^N (y_i-\bar{y})^2}.$ - $ {\frac{{1}}{{2\sigma^2}}}$$ \sum_{{i=1}}^{N}$$ ({y_i-\bar{y}}.$yi - $ \bar{{y}}$$ .{y_i-\bar{y}})^{2}_{}$$ .{-\frac{1}{2\sigma^2}
\sum_{i=1}^N (y_i-\bar{y})^2}]$ = 0

Como o termo exp[...] não pode ser nulo, obtemos
$ {\frac{{d}}{{d\bar{y}}}}$$ [{\sum_{i=1}^N (y_i-\bar{y})^2}.$$ \sum_{{i=1}}^{N}$$ ({y_i-\bar{y}}.$yi - $ \bar{{y}}$$ .{y_i-\bar{y}})^{2}_{}$$ .{\sum_{i=1}^N (y_i-\bar{y})^2}]$ = 0

ou na nossa notação

S = $ \sum_{{i=1}}^{N}$$ ({y_i-\bar{y}}.$yi - $ \bar{{y}}$$ .{y_i-\bar{y}})^{2}_{}$

$ {\frac{dS}{d\bar{y}}=0}$

Continuando com a derivação, obtemos:
0 = $ {\frac{{d}}{{d\bar{y}}}}$$ [{\sum_{i=1}^N (y_i-\bar{y})^2}.$$ \sum_{{i=1}}^{N}$$ ({y_i-\bar{y}}.$yi - $ \bar{{y}}$$ .{y_i-\bar{y}})^{2}_{}$$ .{\sum_{i=1}^N (y_i-\bar{y})^2}]$  
  = $ \sum_{{i=1}}^{N}$ - 2$ ({y_i-\bar{y}}.$yi - $ \bar{{y}}$$ .{y_i-\bar{y}})$  
  = $ \sum_{{i=1}}^{N}$yi - $ \sum_{{i=1}}^{N}$$ \bar{{y}}$  

N$ \bar{{y}}$ = $ \sum_{{i=1}}^{N}$yi

$ \bar{{y}}$ = $ {\frac{{1}}{{N}}}$$ \sum_{{i=1}}^{N}$yi

que é a média simples, como queríamos.

O próximo passo é simplesmente reconhecer que y pode ser uma função, por exemplo:

yi = a + b xi

de modo que

S = $ \sum_{{i=1}}^{N}$$ ({y^o_i - a - b x^o_i}.$yoi - a - b xoi$ .{y^o_i - a - b x^o_i})^{2}_{}$

Minimizando S em relação a a e b, obtemos:

$ {\frac{{dS}}{{da}}}$ = $ \sum_{{i=1}}^{N}$ - 2$ ({y^o_i - a - b\,x^o_i}.$yoi - a - b xoi$ .{y^o_i - a - b\,x^o_i})$ = 0

$ {\frac{{dS}}{{db}}}$ = $ \sum_{{i=1}}^{N}$ - 2xio$ ({y^o_i - a - b\,x^o_i}.$yoi - a - b xoi$ .{y^o_i - a - b\,x^o_i})$ = 0

ou as duas condições:

$ \sum_{{i=1}}^{N}$yio - Na - b$ \sum_{{i=1}}^{N}$xio = 0

$ \sum_{{i=1}}^{N}$xioyio - a$ \sum_{{i=1}}^{N}$xio - b$ \sum_{{i=1}}^{N}$$ ({x_i^o}.$xio$ .{x_i^o})^{2}_{}$ = 0

Em notação matricial podemos escrever as duas condições como:

$ ({\begin{matrix}N&\sum_{i=1}^N x_i^o\cr
\sum_{i=1}^N x_i^o&\sum_{i=1}^N (x_i^o)^2\end{matrix}}.$$ \begin{matrix}N&\sum_{i=1}^N x_i^o\cr
\sum_{i=1}^N x_i^o&\sum_{i=1}^N (x_i^o)^2\end{matrix}$$ .{\begin{matrix}N&\sum_{i=1}^N x_i^o\cr
\sum_{i=1}^N x_i^o&\sum_{i=1}^N (x_i^o)^2\end{matrix}})$$ ({\begin{matrix}a\cr b\cr\end{matrix}}.$$ \begin{matrix}a\cr b\cr\end{matrix}$$ .{\begin{matrix}a\cr b\cr\end{matrix}})$ = $ ({\begin{matrix}\sum_{i=1}^N y_i^o\cr
\sum_{i=1}^N x_i^oy_i^o\cr\end{matrix}}.$$ \begin{matrix}\sum_{i=1}^N y_i^o\cr
\sum_{i=1}^N x_i^oy_i^o\cr\end{matrix}$$ .{\begin{matrix}\sum_{i=1}^N y_i^o\cr
\sum_{i=1}^N x_i^oy_i^o\cr\end{matrix}})$

Para simplificar a notação, vamos definir os colchetes:

$ \sum_{{i=1}}^{N}$xi $ \equiv$ [x]

$ \sum_{{i=1}}^{N}$1 = N $ \equiv$ [1]

desta forma, nossa equação matricial se escreve como:

$ ({\begin{matrix}[1]&[x]\cr [x]&[x^2]\end{matrix}}.$$ \begin{matrix}[1]&[x]\cr [x]&[x^2]\end{matrix}$$ .{\begin{matrix}[1]&[x]\cr [x]&[x^2]\end{matrix}})$$ ({\begin{matrix}a\cr b\end{matrix}}.$$ \begin{matrix}a\cr b\end{matrix}$$ .{\begin{matrix}a\cr b\end{matrix}})$ = $ ({\begin{matrix}[y]\cr [xy]\end{matrix}}.$$ \begin{matrix}[y]\cr [xy]\end{matrix}$$ .{\begin{matrix}[y]\cr [xy]\end{matrix}})$
Estas equações são chamadas de equações normais. Elas podem ser resolvidas com a matriz inversa:
$\left(\begin{matrix}a\cr b\end{matrix}\right)=
\left(\begin{matrix}[1]&[x]\cr [x]&[x^2]\end{matrix}\right)^{-1}
\left(\begin{matrix}[y]\cr [xy]\end{matrix}\right)$
Para sabermos se, quando passamos de uma fitagem de uma reta para uma fitagem de uma parábola, a redução no resíduo, com qualquer normalização, $ S\equiv \chi^2\equiv \sigma^2$ é suficiente para que o termo quadrático seja significativo, podemos definir um parâmetro
$\lambda = \frac{\sigma_{reta}^2 - \sigma_{parab}^2}
{\sigma_{parab}^2}(N-3)$
e determinar o nível de confiabilidade que podemos descartar a hipótese do termo quadrático ser nulo por
$\lambda = F_p(1,n-3)$
onde a distribuição da variável F [Sir Ronald Aylmer Fisher (1890-1962) 1925, Statistical Methods for Research Workers. Oliver and Boyd (Edinburgh)] precisa ser calculada com
$F(a,b)=\frac{\chi^2(a)/a}{\chi^2(b)/b}$
F(a,b) é a variável x da função beta incompleta betai(x,df1,df2), a probabilidade de que uma variável randômica de uma distribuição beta de graus de liberdade df1 e df2 terá valor menor ou igual a x [Statistical Theory and Methodology in Science and Engineering. K. A. Brownlee. London and New York: John Wiley and Sons. 1960].

exemplos de χ2]. Pelo teorema do limite central, se o número de graus de liberdade for grande (maior que 30), a distribuição de χ2 tende à distribuição normal, f=σ1222.

Atente que a definição acima é Fp[(df1-df2),df2], e não F(df1,df2) como usado em geral. Aqui se usou Fp(1,n-3), ou seja df1-df2=1, a diferença no número de graus de liberdade quando passamos de reta para parábola. A definição geral, usada para calcular a probabilidade é

F = [σ2122 * (N-k1)/(N-k2) - 1] * (N-k2)

Outra maneira de analisar a significância do número de parâmetros ajustados é notar que pela definição da distribuição de probabilidades de χ2, o valor esperado de χ2 com k graus de liberdade deve estar no intervalo

⟨χ2⟩±σχ2=k±√(2k)

e k=n-m-1, onde n é o número de observações, m o número de parâmetros ajustados, de modo que

χ2=∑in wi[yi-f(xi,a0,...,am)]2i

tem que estar neste intervalo. Se for menor, o número de parâmetros está superestimado.

Mínimos quadrados lineares

Os mínimos quadrados são lineares quando podemos resolver as equações normais usando álgebra linear. O exemplo dado acima é um mínimo quadrado linear, pois podemos resolver as equações para a e b exatamente.

Mínimos quadrados não lineares

Muitos problemas interessantes não podem ser resolvidos linearmente. Por exemplo, qual é o melhor valor de $ \alpha$ em

y = exp(- $ \alpha$x)

Pela nossa definição de S:

S=$ \sum_{{i=1}}^{N}$$ [{y_i^o - \exp (-\alpha x_i)}.$yio - exp$ ({-\alpha x_i}.$ - $ \alpha$xi$ .{-\alpha x_i})$$ .{y_i^o - \exp (-\alpha x_i)}]^{2}_{}$

e quando minimizamos S:

0 = $ {\frac{{dS}}{{d\alpha}}}$ = $ \sum_{{i=1}}^{N}$2$ [{y_i^o - \exp (-\alpha x_i)}.$yio - exp$ ({-\alpha x_i}.$ - $ \alpha$xi$ .{-\alpha x_i})$$ .{y_i^o - \exp (-\alpha x_i)}]$$ [{- \exp (-\alpha x_i)}.$ -exp$ ({-\alpha x_i}.$ - $ \alpha$xi$ .{-\alpha x_i})$$ .{- \exp (-\alpha x_i)}]$$ ({-x_i}.$ - xi$ .{-x_i})$

ou seja

0 = $ \sum_{{i=1}}^{N}$yioxiexp$ ({-\alpha x_i}.$ - $ \alpha$xi$ .{-\alpha x_i})$ - $ \sum_{{i=1}}^{N}$xiexp$ ({-2\alpha x_i}.$ -2$ \alpha$xi$ .{-2\alpha x_i})$

Que podemos escrever, na notação de colchetes, como:

0 = $ [{xy\exp(-\alpha x)}.$xy exp$ ({-\alpha x}.$ - $ \alpha$x$ .{-\alpha x})$$ .{xy\exp(-\alpha x)}]$ - $ [{x\exp(-2\alpha x)}.$x exp$ ({-2\alpha x}.$ -2$ \alpha$x$ .{-2\alpha x})$$ .{x\exp(-2\alpha x)}]$

Esta equação não pode ser resolvida usando-se álgebra linear.

Precisamos utilizar técnicas diferentes. A técnica mais empregada e, de fato, a técnica que se utiliza quando chamamos de mínimos quadrados não lineares, é a linearização do problema.

A idéia, aplicada ao problema acima, é:

y = exp(- $ \alpha$x)

Escolhemos um valor inicial de $ \alpha$, chamado $ \alpha_{0}^{}$, e definimos:

$ \alpha$ = $ \alpha_{0}^{}$ + $ \Delta$$ \alpha$

Definindo

y0 = exp(- $ \alpha_{0}^{}$x)

Em primeira ordem, isto é, linear:

y = y0 + $ {\frac{{dy}}{{d\alpha}}}$$ \vert _{{\alpha_0}}^{}$$ \Delta$$ \alpha$

y = exp$ ({-\alpha_0 x}.$ - $ \alpha_{0}^{}$x$ .{-\alpha_0 x})$ - x exp$ ({-\alpha_0 x}.$ - $ \alpha_{0}^{}$x$ .{-\alpha_0 x})$$ \Delta$$ \alpha$

Agora y é linear em $ \Delta$$ \alpha$ e podemos usar os mínimos quadrados lineares para encontrar a correção $ \Delta$$ \alpha$:

S $ \equiv$ $ \sum_{{i=1}}^{N}$$ ({y^o_i - y_i}.$yoi - yi$ .{y^o_i - y_i})^{2}_{}$

que se torna

S $ \equiv$ $ \sum_{{i=1}}^{N}$$ ({y^o_i - y_{0,i}-\frac{dy}{d\alpha}\vert _{\alpha_0}\Delta \alpha}.$yoi - y0, i - $ {\frac{{dy}}{{d\alpha}}}$$ \vert _{{\alpha_0}}^{}$$ \Delta$$ \alpha$$ .{y^o_i - y_{0,i}-\frac{dy}{d\alpha}\vert _{\alpha_0}\Delta \alpha})^{2}_{}$

que minimizando

0 = $ {\frac{{dS}}{{d(\Delta\alpha)}}}$ = $ \sum_{{i=1}}^{N}$ - 2$ {\frac{{dy}}{{d\alpha}}}$$ \vert _{{\alpha_0}}^{}$$ ({y^o_i - y_{0,i}-\frac{dy}{d\alpha}\vert _{\alpha_0}\Delta \alpha}.$yoi - y0, i - $ {\frac{{dy}}{{d\alpha}}}$$ \vert _{{\alpha_0}}^{}$$ \Delta$$ \alpha$$ {y^o_i - y_{0,i}-\frac{dy}{d\alpha}\vert _{\alpha_0}
\Delta \alpha})$

0 = $ \sum_{{i=1}}^{N}$$ [{y^o_i \frac{dy}{d\alpha}\vert _{\alpha_0}
-(\frac{dy}{d\alpha}\vert _{\alpha_0})^2\Delta \alpha}.$yoi$ {\frac{{dy}}{{d\alpha}}}$$ \vert _{{\alpha_0}}^{}$ - y0, i$ {\frac{{dy}}{{d\alpha}}}$$ \vert _{{\alpha_0}}^{}$ - $ ({\frac{dy}{d\alpha}\vert _{\alpha_0}}.$$ {\frac{{dy}}{{d\alpha}}}$$ \vert _{{\alpha_0}}^{}$$ .{\frac{dy}{d\alpha}\vert _{\alpha_0}})^{2}_{}$$ \Delta$$ \alpha$$ .{y^o_i \frac{dy}{d\alpha}\vert _{\alpha_0}
-y_...
...0}
-(\frac{dy}{d\alpha}\vert _{\alpha_0})^2\Delta \alpha}]$

ou, na notação dos colchetes:
0 = $ [{y^o \frac{dy}{d\alpha}}.$yo$ {\frac{{dy}}{{d\alpha}}}$$ .{y^o \frac{dy}{d\alpha}}]$ - $ [{y_0 \frac{dy}{d\alpha}}.$y0$ {\frac{{dy}}{{d\alpha}}}$$ .{y_0 \frac{dy}{d\alpha}}]$ - $ [{(\frac{dy}{d\alpha})^2\Delta\alpha}.$$ ({\frac{dy}{d\alpha}}.$$ {\frac{{dy}}{{d\alpha}}}$$ .{\frac{dy}{d\alpha}})^{2}_{}$$ \Delta$$ \alpha$$ .{(\frac{dy}{d\alpha})^2\Delta\alpha}]$

Que podemos resolver para $ \Delta$$ \alpha$:
$ \Delta$$ \alpha$ = $ {\frac{{[y^o \frac{dy}{d\alpha}]
-[y_0 \frac{dy}{d\alpha}]}}{{[(\frac{dy}{d\alpha})^2]}}}$

e finalmente obter o valor revisado de $ \alpha$:
$ \alpha$ = $ \alpha_{0}^{}$ + $ \Delta$$ \alpha$

Note entretanto que o valor revisado de $ \alpha$ não é o melhor valor de $ \alpha$, mas somente uma melhor aproximação do que $ \alpha_{0}^{}$. Isto ocorre porque $ \Delta$$ \alpha$ é a solução do problema linearizado e não do problema real. Portanto, precisamos iterar, isto é, utilizar este valor revisado de $ \alpha$ como um novo $ \alpha_{0}^{}$ e obter um novo valor revisado.

Formulação Geral

Se a função y for uma função de k parâmetros que queremos ajustar:
yi = y$ ({x_i,a_1,a_2,\ldots,a_k}.$xi, a1, a2,..., ak$ .{x_i,a_1,a_2,\ldots,a_k})$

colocamos
$ \begin{matrix}
a_1=a_{1,0}+\Delta a_1\\
a_2=a_{2,0}+\Delta a_2\\
\vdots\\
a_k=a_{k,0}+\Delta a_k
\end{matrix}$

com a hipótese de que $ \Delta$ai $ \ll$ ai, 0.

Então podemos linearizar

yi = y$ ({x_i,a_{1,0},a_{2,0},\ldots,a_{k,0}}.$xi, a1, 0, a2, 0,..., ak, 0$ .{x_i,a_{1,0},a_{2,0},\ldots,a_{k,0}})$ + $ {\frac{{dy}}{{da_1}}}$$ \vert _{{a_n=a_{n,0}}}^{}$$ \Delta$a1 + $ {\frac{{dy}}{{da_2}}}$$ \vert _{{a_n=a_{n,0}}}^{}$$ \Delta$a2 + ...  
    + ... + $ {\frac{{dy}}{{da_k}}}$$ \vert _{{a_n=a_{n,0}}}^{}$$ \Delta$ak  

notando que as derivadas são calculadas para todos an = an, 0.

Chamando

yi, 0 = y$ ({x_i,a_{1,0},a_{2,0},\ldots,a_{k,0}}.$xi, a1, 0, a2, 0,..., ak, 0$ .{x_i,a_{1,0},a_{2,0},\ldots,a_{k,0}})$
e
$ {\frac{{dy_i}}{{da_j}}}$$ \vert _{{a_0}}^{}$$ \Delta$aj = $ {\frac{{dy_i}}{{da_j}}}$$ \vert _{{a_n=a_{n,0}}}^{}$$ \Delta$aj
podemos escrever
yi = yi, 0 + $ \sum_{{j=1}}^{k}$$ {\frac{{dy_i}}{{da_j}}}$$ \vert _{{a_0}}^{}$$ \Delta$aj

onde o subscrito i significa calculado no ponto xi.

Podemos agora calcular S:

S $ \equiv$ $ \sum_{{i=1}}^{N}$$ ({y^o_i - y_i}.$yoi - yi$ .{y^o_i - y_i})^{2}_{}$
S $ \equiv$ $ \sum_{{i=1}}^{N}$$ ({y^o_i - y_{i,0} -
\sum_{j=1}^k \frac{dy_i}{da_j}\vert _{a_0}\Delta a_j}.$yoi - yi, 0 - $ \sum_{{j=1}}^{k}$$ {\frac{{dy_i}}{{da_j}}}$$ \vert _{{a_0}}^{}$$ \Delta$aj$ .{y^o_i - y_{i,0} -
\sum_{j=1}^k \frac{dy_i}{da_j}\vert _{a_0}\Delta a_j})^{2}_{}$

que minimizando com respeito a $ \Delta$am:
0 = $ {\frac{{dS}}{{d(\Delta a_m)}}}$ = $ \sum_{{i=1}}^{N}$2$ ({y^o_i - y_{i,0} -
\sum_{j=1}^k \frac{dy_i}{da_j}\vert _{a_0}\Delta a_j}.$yoi - yi, 0 - $ \sum_{{j=1}}^{k}$$ {\frac{{dy_i}}{{da_j}}}$$ \vert _{{a_0}}^{}$$ \Delta$aj$ .{y^o_i - y_{i,0} -
\sum_{j=1}^k \frac{dy_i}{da_j}\vert _{a_0}\Delta a_j})$$ ({-\frac{dy_i}{da_m}\vert _{a_0}}.$ - $ {\frac{{dy_i}}{{da_m}}}$$ \vert _{{a_0}}^{}$$ .{-\frac{dy_i}{da_m}\vert _{a_0}})$
0 = $ \sum_{{i=1}}^{N}$$ ({y^o_i - y_{i,0}}.$yoi - yi, 0$ .{y^o_i - y_{i,0}})$$ ({\frac{dy_i}{da_m}\vert _{a_0}}.$$ {\frac{{dy_i}}{{da_m}}}$$ \vert _{{a_0}}^{}$$ .{\frac{dy_i}{da_m}\vert _{a_0}})$ - $ \sum_{{j=1}}^{k}$$ {\frac{{dy_i}}{{da_j}}}$$ \vert _{{a_0}}^{}$$ \Delta$aj$ ({\sum_{i=1}^N\frac{dy_i}{da_m}\vert _{a_0}}.$$ \sum_{{i=1}}^{N}$$ {\frac{{dy_i}}{{da_m}}}$$ \vert _{{a_0}}^{}$$ .{\sum_{i=1}^N\frac{dy_i}{da_m}\vert _{a_0}})$

que na notação dos colchetes pode ser escrita como:
0 = $ [{(y^o - y_0)\frac{dy}{da_m}}.$$ ({y^o - y_0}.$yo - y0$ .{y^o - y_0})$$ {\frac{{dy}}{{da_m}}}$$ .{(y^o - y_0)\frac{dy}{da_m}}]$ - $ \sum_{{j=1}}^{k}$$ [{\frac{dy}{da_m}\frac{dy}{da_j}}.$$ {\frac{{dy}}{{da_m}}}$$ {\frac{{dy}}{{da_j}}}$$ .{\frac{dy}{da_m}\frac{dy}{da_j}}]$$ \Delta$aj

Ou, em notação matricial
$ ({\begin{matrix}
[\frac{dy}{da_1}\frac{dy}{da_...
...}]&\cdots
[\frac{dy}{da_k}\frac{dy}{da_k}]\end{matrix}}.$$ \begin{matrix}
[\frac{dy}{da_1}\frac{dy}{da_1}]&
...
...dy}{da_2}]&\cdots
[\frac{dy}{da_k}\frac{dy}{da_k}]\end{matrix}$$ .{\begin{matrix}
[\frac{dy}{da_1}\frac{dy}{da_...
...}]&\cdots
[\frac{dy}{da_k}\frac{dy}{da_k}]\end{matrix}})$$ ({\begin{matrix}
\Delta a_1\\
\Delta a_2\\
\vdots\\
\Delta a_k\end{matrix}}.$$ \begin{matrix}
\Delta a_1\\
\Delta a_2\\
\vdots\\
\Delta a_k\end{matrix}$$ .{\begin{matrix}
\Delta a_1\\
\Delta a_2\\
\vdots\\
\Delta a_k\end{matrix}})$ = $ ({\begin{matrix}
[(y^o - y_0)\frac{...
...ots\\
[(y^o - y_0)\frac{dy}{da_k}]
\end{matrix}}.$$ \begin{matrix}
[(y^o - y_0)\frac{dy}{da_1}]...
...\\
\vdots\\
[(y^o - y_0)\frac{dy}{da_k}]
\end{matrix}$$ .{\begin{matrix}
[(y^o - y_0)\frac{...
...ots\\
[(y^o - y_0)\frac{dy}{da_k}]
\end{matrix}})$

Esta equação matricial pode agora ser revolvida por álgebra matricial para encontrar as correções $ \Delta$am.

No caso linear

$\left(\begin{matrix}
\left[\frac{dy}{da_1}\frac{dy}{da_1}\right] ...\\
\left[\left(y^o - y_0\right)\frac{dy}{da_k}\right]
\end{matrix}\right)$

Determinação das incertezas e Propagação de Erros

Somente um valor com sua incerteza permite reconstruir a distribuição de probabilidades de um parâmetro. A maneira correta de determinar as incertezas nos parâmetros é calculando a variança de um parâmetro qualquer z:
$ \sigma_{z}^{2}$ $ \equiv$ $ \langle$ ($ \Delta$z)2 $ \rangle$
onde o quadrado é necessário pelas mesmas considerações do valor de S, $ \langle$a $ \rangle$ significa a média e onde
$ \Delta$z = zcalculado - zobservado
Agora suponhamos que z seja uma função de duas variáveis x e y, z = z(x, y). Podemos expandir por série de Taylor [Brook Taylor (1685-1731), Methodus incrementorum directa et inversa (1715)]:
$ \Delta$z = $ {\frac{{dz}}{{dx}}}$$ \Delta$x + $ {\frac{{dz}}{{dy}}}$$ \Delta$y
de onde obtemos:
$ \sigma_{z}^{2}$ $ \equiv$ $ \langle$($ \Delta$z)2$ \rangle$ = $ g\langle$$ ({\frac{dz}{dx}\Delta x + \frac{dz}{dy}\Delta y}.$$ {\frac{{dz}}{{dx}}}$$ \Delta$x + $ {\frac{{dz}}{{dy}}}$$ \Delta$y$ .{\frac{dz}{dx}\Delta x + \frac{dz}{dy}\Delta y})^{2}_{}$$ g\rangle$  
  = $ g\langle$$ ({\frac{dz}{dx}}.$$ {\frac{{dz}}{{dx}}}$$ .{\frac{dz}{dx}})^{2}_{}$$ ({\Delta x}.$$ \Delta$x$ .{\Delta x})^{2}_{}$ + 2$ {\frac{{dz}}{{dx}}}$$ {\frac{{dz}}{{dy}}}$$ \Delta$x$ \Delta$y + $ ({\frac{dz}{dy}}.$$ {\frac{{dz}}{{dy}}}$$ .{\frac{dz}{dy}})^{2}_{}$$ ({\Delta y}.$$ \Delta$y$ .{\Delta y})^{2}_{}$$ g\rangle$  
Se as variáveis x e y são separáveis, podemos reduzir a equação acima a
$ \sigma_{z}^{2}$ = $ ({\frac{dz}{dx}}.$$ {\frac{{dz}}{{dx}}}$$ .{\frac{dz}{dx}})^{2}_{}$$ \big\langle$$ ({\Delta x}.$$ \Delta$x$ .{\Delta x})^{2}_{}$$ \big\rangle$ + 2$ {\frac{{dz}}{{dx}}}$$ {\frac{{dz}}{{dy}}}$$ \big\langle$$ \Delta$x$ \Delta$y$ \big\rangle$ + $ ({\frac{dz}{dy}}.$$ {\frac{{dz}}{{dy}}}$$ .{\frac{dz}{dy}})^{2}_{}$$ \big\langle$$ ({\Delta y}.$$ \Delta$y$ .{\Delta y})^{2}_{}$$ \big\rangle$
e, por definição:
$ \sigma_{x}^{2}$ $ \equiv$ $ \langle$($ \Delta$x)2$ \rangle$
$ \sigma_{y}^{2}$ $ \equiv$ $ \langle$($ \Delta$y)2$ \rangle$
$ \sigma_{{xy}}^{2}$ $ \equiv$ $ \langle$$ \Delta$x$ \Delta$y$ \rangle$
de modo que
$ \sigma_{z}^{2}$ = $ ({\frac{dz}{dx}}.$$ {\frac{{dz}}{{dx}}}$$ .{\frac{dz}{dx}})^{2}_{}$$ \sigma_{x}^{2}$ + 2$ {\frac{{dz}}{{dx}}}$$ {\frac{{dz}}{{dy}}}$$ \sigma_{{xy}}^{2}$ + $ ({\frac{dz}{dy}}.$$ {\frac{{dz}}{{dy}}}$$ .{\frac{dz}{dy}})^{2}_{}$$ \sigma_{y}^{2}$

E como obtemos $ \sigma_{x}^{}$, $ \sigma_{y}^{}$ e $ \sigma_{{xy}}^{}$? Definindo a matriz covariança.

Matriz Covariança

Definimos a matriz covariança como

COV = $ ({\begin{matrix}\sigma_x^2&\sigma_{xy}^2\\  \sigma_{xy}^2&\sigma_y^2\end{matrix}}.$$ \begin{matrix}\sigma_x^2&\sigma_{xy}^2\\  \sigma_{xy}^2&\sigma_y^2\end{matrix}$$ .{\begin{matrix}\sigma_x^2&\sigma_{xy}^2\\  \sigma_{xy}^2&\sigma_y^2\end{matrix}})$

De modo que, se a equação normal matricial que usamos para achar os coeficientes é escrita como

$ \bf {MA=Y}$

de modo que a solução é dada por

$ \bf A = M^{-1}Y$

onde $ \bf {M^{-1}}$ é a matriz inversa da matriz $ \bf M$, vemos que a matriz covariança é dada por
$ \bf {COV(A) = \sigma^2 M^{-1}}$
onde
$ \sigma^{2}_{}$ = $ {\frac{{S}}{{N-k}}}$ = $ {\frac{{\sum_{i=1}^N (y^o_i - y_i)^2}}{{N-k}}}$

Desta maneira é muito fácil calcular as incertezas nos parâmetros, pois somente precisamos multiplar a matriz inversa que usamos para obter os valores dos parâmetros, por $ \sigma^{2}_{}$, obtendo a matriz covariança.

χ2

Mas o que fazemos se a função a ser fitada não é analítica, como por exemplo, um espectro resultante de um modelo de atmosferas? Er-Ho Zhang, Edward L. Robinson e R. Edward Nather (1986, Astrophysical Journal, 305, 740) demonstraram que no caso medirmos vários parâmetros simultaneamente, onde a mudança de alguns pode gerar mudanças significativas, enquanto que a mudança de outros parâmetros gera pequenas diferenças e, ainda, que os parâmetros não são independentes, isto é, estão correlacionados, precisamos maximizar a probabilidade (likelihood) de que os parâmetros x1,..., xk sejam medidos pelos valores x1o,..., xko, definindo incertezas $ \epsilon_{i}^{}$ = xio - xi,
f ($ \epsilon_{1}^{}$,...,$ \epsilon_{k}^{}$) = $ {\frac{{1}}{{(2\pi)^{k/2}\sigma_1\cdots \sigma_k}}}$exp$ [{-\frac{1}{2}(\frac{\epsilon_1^2}{\sigma_1^2}
+\cdots + \frac{\epsilon_k^2}{\sigma_k^2})}.$ - $ {\frac{{1}}{{2}}}$$ ({\frac{\epsilon_1^2}{\sigma_1^2}
+\cdots + \frac{\epsilon_k^2}{\sigma_k^2}}.$$ {\frac{{\epsilon_1^2}}{{\sigma_1^2}}}$ + ... + $ {\frac{{\epsilon_k^2}}{{\sigma_k^2}}}$$ .{\frac{\epsilon_1^2}{\sigma_1^2}
+\cdots + \frac{\epsilon_k^2}{\sigma_k^2}})$$ .{-\frac{1}{2}(\frac{\epsilon_1^2}{\sigma_1^2}
+\cdots + \frac{\epsilon_k^2}{\sigma_k^2})}]$
e, portanto, maximizarmos a probabilidade é equivalente a minimizarmos
$\chi^2 \equiv \tilde{S} \equiv \frac{\epsilon_1^2}{\sigma_1^2}
+\cdots + \frac{\epsilon_k^2}{\sigma_k^2}$
Note que o valor de $ \tilde{{S}}$ tem normalização diferente da de S definido na equação (1). A normalização diferente é que distingue χ2 de S2.

Se medimos os valores Ioi, por exemplo, o espectro, e calculamos Ii com os parâmetros x1,..., xk, assumindo que a distribuição de erros é normal e que todas as varianças são iguais ($ \sigma_{1}^{}$ = $ \sigma_{2}^{}$ = ... = $ \sigma_{k}^{}$ = W), obtemos

$\tilde{S} = \frac{\sum_i^N \left(I_i^o - I_i\right)^2}{W^2}$
de modo que, se $ \tilde{{S}}_{0}^{}$ é o valor mínimo de $ \tilde{{S}}$, e se mudarmos o valor do parâmetro i por um delta d, mantendo todos os outros parâmetros constantes, obteremos um novo valor de $ \tilde{{S}}$
$ \tilde{{S}}$ = $ \tilde{{S}}_{0}^{}$ + d2/$ \sigma_{{ii}}^{}$
e, portanto,
$ \sigma_{i}^{2}$ = $ \sigma_{{ii}}^{}$ = $ {\frac{{d^2}}{{\tilde{S}-\tilde{S}_0}}}$ (2)
Portanto podemos obter a incerteza em cada parâmetro minimizando $ \tilde{{S}}$ com todos os parâmetros livres, para obter $ \tilde{{S}}_{0}^{}$, fixando o valor do parâmetro xi para o valor que minimizou $ \tilde{{S}}$ mais um delta, digamos 5%, e minimizando novamente $ \tilde{{S}}$ com este valor de xi fixo, obtendo um novo $ \tilde{{S}}$. Com a equação (2) podemos então estimar o valor de sua incerteza, e repetir o processo para os outros parâmetros.

Podemos encontrar o termo de correlação usando uma transformação de variáveis

$d_i\prime = (d_i+d_j)/\sqrt{2}$
$d_j\prime = (d_i-d_j)/\sqrt{2}$
com $ d_i=d_j=d=d_{i\prime}/\sqrt{2}$. Pela propagação de erros,
$\sigma_{i\prime}^2=\frac{1}{2}(\sigma_i^2 + 2\sigma_{ij} + \sigma_j^2)$
de modo que
$\sigma_{ij} = \sigma_{i\prime}^2-\frac{1}{2}(\sigma_i^2 + \sigma_j^2) = \frac{d^2}{\tilde{S}-\tilde{S}_0}-\frac{1}{2}(\sigma_i^2 + \sigma_j^2)$
Se a normalização não for unitária, devemos usar:
sigma_i^2 = \sigma_{ii} = \frac{S}{n-k+1}\frac{d^2}{\tilde{S}-\tilde{S}_0}$ (1)
onde n é o número de pontos, k o número de parâmetros fitados, e o +1 é porque um dos parâmetros foi mantido fixo, reduzindo o número de graus de liberdade.

Como demonstrado por Edward Lewis Robinson (1945-), no seu livro de 2016, Data Analysis for Scientists and Engineers, Princeton University Press, a distribuição de probabilidades de χ2 para k graus de liberdade é dada por

$f_k(\chi^2) = \frac{(\chi^2)^{(k-2)/2}}{2^{k/2}(k/2-1)!}
\exp{[-\frac{1}{2}\chi^2]}$
distribuicao de chi2
Note que o valor médio de χ2 para k graus de liberdade é k, a moda é k-2, e a variança σ2χ2=2k.
Algumas vezes se usa um χ2 reduzido de modo que sua média seja 1:
$\chi^2_{red} \equiv \frac{1}{k}\chi^2 = \frac{1}{k}\sum_i^k
\frac{\epsilon_i^2}{\sigma_i^2}$
Embora o valor médio do χ2red seja 1, a variança continua dependente do número de graus de liberdade k
⟨χ2red⟩=1
σ2χ2red=2/k
Sabemos que o teorema do limite central garante que a distribuição de χ2 para k graus de liberdade se aproxima da distribuição gaussiana para grandes valores de k. Esta aproximação é boa para k>30. O grande problema do método χ2 é que ele assume que as incertezas das medidas são bem determinadas e não correlacionadas. Se as incertezas forem subestimadas, o valor de χ2 não tem significado.

Estimativa Robusta

O grande problema do método de mínimos quadrados é a grande influência de pontos com resíduo muito grande, os outliers. Para se reduzir a influência destes resíduos muito grandes, pode-se minimizar funções que não crescem tão rapidamente quanto S2, como por exemplo, mimizando a função discrepância, definida como:

f(S)=\left\{
\begin{array}{ll}
S^2&\mbox{se $\vert S\vert\leq c$}\\
c^2&\mbox{se $\vert S\vert> c$}
\end{array}
onde c é uma constante.

Probabilidade

A função normal, centrada em x=0 e com largura em 1/e em x=1 e x=-1, é definida como

P(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^2}

Probabilidade
e sua integral é a distribuição normal
f(x)=\frac{1}{\sqrt{2\pi}}\int^x_{-\infty} e^{-\frac{1}{2}y^2} dy

Com esta definição
\int^{+\infty}_{-\infty}P(x)dx=1
Se tivermos uma distribuição centrada em $x_o$, e de largura $\sigma^2$,

P(x-x_o)=\frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x-x_o)^2}{2\sigma^2}}
O valor da integral de x até infinito, centrada em 0 e com desvio padrão 1, é
x=0P(>x)=0,5
x=2P(>x)=0,02275
x=3P(>x)=0,00135
x=4P(>x)=0,00003167
x=5P(>x)=2,867 × 10-7
x=6P(>x)=9,866 × 10-10
x=7P(>x)=1,280 × 10-12
x=8P(>x)=6,221 × 10-16
x=9P(>x)=1,129 × 10-19
Se a distribuição é discreta, com N valores:
\sum_{i=0}^N P(x_i) \equiv 1

de modo que precisamos normalizar nossas probabilidades
P(x-x_o)=\frac{1}{\sum_{i=0}^N P(x_i)}\frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{(x-x_o)^2}{2\sigma^2}}

Distribuição Binominal

Seja um problema que só admite duas soluções, como ao jogarmos uma moeda. Chamemos de b(k;n,p) a probabilidade que n tentativas com probabilidade intrínseca de sucesso p, e q=1-p a probabilidade intrínseca de insucesso, resulte em k sucessos e n-k insucessos. Esta probabilidade é dada pela distribuição binominal
b(k;n,p) = {n\choose k} p^k q^{n-k}=\frac{1}{\sqrt{npq}}P\(\frac{k-np}{\sqrt{npq}})
Gaussiana 2D
Distribuição gaussiana (normal) para duas variáveis não correlacionadas.
Para duas variáveis correlacionadas, x1 medida a partir de μ1 e x2 medida a partir de μ2:

P(x_1,x_2)=\frac{1}{2\pi\sigma_1\sigma_2\sqrt{(1-\rho^2)}}\exp{[-\frac{z}{2(1-\rho^2)}]}

onde

z=\frac{(x_1-\mu_1)^2}{\sigma_1^2} - 2\rho\frac{(x_1-\mu_1)(x_2-\mu_2)}{\sigma_1\sigma_2} + \frac{(x_2-\mu_2)^2}{\sigma_2^2}

\rho=cor(x_1,x_2)=\frac{\sigma_{1,2}}{\sigma_1\sigma_2}


Variança da média
Pesos nas medidas
Cálculo de Distribuição Normal
Outras Distribuições
Data analysis recipes: Fitting a model to data, David W. Hogg, Jo Bovy (NYU), Dustin Lang, 2010. arXiv:1008.4686 Tutorial sobre diferentes métodos de ajuste linear a dados, com os problemas de excluir outliers modificando a probabilidade, explicitando que o método de mínimos quadrados exclui incertezas na variável de medida (x), isto é as incertezas são somente em uma dimensão, e assume que as incertezas são todas gaussianas. Inclui exercícios e uma introdução à estatítica Baesiana. Fornece uma visão crítica sobre os métodos comumente adotados. É interessante por exemplo a crítica dos autores à metodologia usada para derivar a relação de Tully-Fisher, e a crítica ao uso de PCA (principal component analysis), em ambos casos por terem incertezas nas duas dimensões.

Error estimation in astronomy: A guide. Rene Andrae. 2010. arXiv:1009.2755 Tutorial sobre estimativa das incertezas em parâmetros.

Dos and don'ts of reduced chi-squared (χ2), Rene Andrae, Tim Schulze-Hartung, Peter Melchior, arXiv:1012.3754.

No Excel, mínimos quadrados lineares são fitados em Tools-Data Analysis-Regression. Se este pacote não estiver instalado, ele pode ser adicionado com Tools-Add Ins-Analysis Tool Pack e Solver. O OpenOffice Calc também pode ser usado. Com uma subrotina em Fortran de inversão de matrizes, se pode facilmente implementar um algoritmo de mínimos quadrados.

Maiores detalhes e todas as derivações estão no livro Data Analysis for Scientists and Engineers, de Edward L. Robinson, 2016, Princeton University Press. Statistical Calculations
Handbook of Bio Statistics
Fitando períodos com Lorentzianas - Keaton
Relação entre amplitude média e sigma em uma transformada de Fourier

©
Modificada em 20 maio 2021