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

y^o_i = valores observados de y

y_i = 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.

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).

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 {y_i}. A probabilidade de observar este conjunto é dada por

$$({y_1,y_2,y_3,\ldots,y_N}. .{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}}}} [{-\frac{1}{2\sigma^2} (y_1-\bar{y})^2}. {\frac{{1}}{{2\sigma^2}}} ({y_1-\bar{y}}. \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}}}} [{-\frac{1}{2\sigma^2} (y_2-\bar{y})^2}. {\frac{{1}}{{2\sigma^2}}} ({y_2-\bar{y}}. \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}}}} [{-\frac{1}{2\sigma^2} (y_N-\bar{y})^2} {\frac{{1}}{{2\sigma^2}}} ({y_N-\bar{y}}. \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}_{} [{-\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}}. \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 {y_i}, 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}.$y₁, y₂, y₃,..., y_N$.{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}}.$y_i - $\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}}.$y_i - $\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}}.$y_i - $\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}}.$y_i - $\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}}.$y_i - $\bar{{y}}$$.{y_i-\bar{y}})^{2}_{}$

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

Continuando com a derivação, obtemos:

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

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

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

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

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

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

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

y_i = a + b x_i

de modo que

S = $\sum_{{i=1}}^{N}$$({y^o_i - a - b x^o_i}.$y^o_i - a - b x^o_i$.{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}.$y^o_i - a - b x^o_i$.{y^o_i - a - b\,x^o_i})$ = 0

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

ou as duas condições:

$\sum_{{i=1}}^{N}$y_i^o - Na - b$\sum_{{i=1}}^{N}$x_i^o = 0

$\sum_{{i=1}}^{N}$x_i^oy_i^o - a$\sum_{{i=1}}^{N}$x_i^o - b$\sum_{{i=1}}^{N}$$({x_i^o}.$x_i^o$.{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}$x_i $\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 χ²]. Pelo teorema do limite central, se o número de graus de liberdade for grande (maior que 30), a distribuição de χ² tende à distribuição normal, f=σ₁²/σ₂².

Atente que a definição acima é F_p[(df1-df2),df2], e não F(df1,df2) como usado em geral. Aqui se usou F_p(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 = [σ²₁/σ²₂ * (N-k₁)/(N-k₂) - 1] * (N-k₂)

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 χ², o valor esperado de χ² com k graus de liberdade deve estar no intervalo

⟨χ²⟩±σ_χ²=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

χ²=∑_iⁿ w_i[y_i-f(x_i,a₀,...,a_m)]²/σ_i

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)}.$y_i^o - exp$({-\alpha x_i}.$ - $\alpha$x_i$.{-\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)}.$y_i^o - exp$({-\alpha x_i}.$ - $\alpha$x_i$.{-\alpha x_i})$$.{y_i^o - \exp (-\alpha x_i)}]$$[{- \exp (-\alpha x_i)}.$ -exp$({-\alpha x_i}.$ - $\alpha$x_i$.{-\alpha x_i})$$.{- \exp (-\alpha x_i)}]$$({-x_i}.$ - x_i$.{-x_i})$

ou seja

0 = $\sum_{{i=1}}^{N}$y_i^ox_iexp$({-\alpha x_i}.$ - $\alpha$x_i$.{-\alpha x_i})$ - $\sum_{{i=1}}^{N}$x_iexp$({-2\alpha x_i}.$ -2$\alpha$x_i$.{-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

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

Em primeira ordem, isto é, linear:

y = y₀ + ${\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}.$y^o_i - y_i$.{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}.$y^o_i - y_{0, 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}.$y^o_i - y_{0, 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}.$y^o_i${\frac{{dy}}{{d\alpha}}}$$\vert _{{\alpha_0}}^{}$ - y_{0, 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}}.$y^o${\frac{{dy}}{{d\alpha}}}$$.{y^o \frac{dy}{d\alpha}}]$ - $[{y_0 \frac{dy}{d\alpha}}.$y₀${\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:

y_i = y$({x_i,a_1,a_2,\ldots,a_k}.$x_i, a₁, a₂,..., a_k$.{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$a_i $\ll$ a_{i, 0}.

Então podemos linearizar

$$({x_i,a_{1,0},a_{2,0},\ldots,a_{k,0}}. .{x_i,a_{1,0},a_{2,0},\ldots,a_{k,0}}) {\frac{{dy}}{{da_1}}} \vert _{{a_n=a_{n,0}}}^{} \Delta {\frac{{dy}}{{da_2}}} \vert _{{a_n=a_{n,0}}}^{} \Delta$$ y_i =

$${\frac{{dy}}{{da_k}}} \vert _{{a_n=a_{n,0}}}^{} \Delta$$

notando que as derivadas são calculadas para todos a_n = a_{n, 0}.

Chamando

y_{i, 0} = y$({x_i,a_{1,0},a_{2,0},\ldots,a_{k,0}}.$x_i, a_{1, 0}, a_{2, 0},..., a_{k, 0}$.{x_i,a_{1,0},a_{2,0},\ldots,a_{k,0}})$

e

${\frac{{dy_i}}{{da_j}}}$$\vert _{{a_0}}^{}$$\Delta$a_j = ${\frac{{dy_i}}{{da_j}}}$$\vert _{{a_n=a_{n,0}}}^{}$$\Delta$a_j

podemos escrever

y_i = y_{i, 0} + $\sum_{{j=1}}^{k}$${\frac{{dy_i}}{{da_j}}}$$\vert _{{a_0}}^{}$$\Delta$a_j

onde o subscrito i significa calculado no ponto x_i.

Podemos agora calcular S:

S $\equiv$ $\sum_{{i=1}}^{N}$$({y^o_i - y_i}.$y^o_i - y_i$.{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}.$y^o_i - y_{i, 0} - $\sum_{{j=1}}^{k}$${\frac{{dy_i}}{{da_j}}}$$\vert _{{a_0}}^{}$$\Delta$a_j$.{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$a_m:

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}.$y^o_i - y_{i, 0} - $\sum_{{j=1}}^{k}$${\frac{{dy_i}}{{da_j}}}$$\vert _{{a_0}}^{}$$\Delta$a_j$.{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}}.$y^o_i - y_{i, 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$a_j$({\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}.$y^o - y₀$.{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$a_j

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

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)² $\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 = z_calculado - z_observado

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 \rangle g\langle ({\frac{dz}{dx}\Delta x + \frac{dz}{dy}\Delta y}. {\frac{{dz}}{{dx}}} \Delta {\frac{{dz}}{{dy}}} \Delta .{\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 .{\Delta x})^{2}_{} {\frac{{dz}}{{dx}}} {\frac{{dz}}{{dy}}} \Delta \Delta ({\frac{dz}{dy}}. {\frac{{dz}}{{dy}}} .{\frac{dz}{dy}})^{2}_{} ({\Delta y}. \Delta .{\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)²$\rangle$

$\sigma_{y}^{2}$ $\equiv$ $\langle$($\Delta$y)²$\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.

χ²

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 x₁,..., x_k sejam medidos pelos valores x₁^o,..., x_k^o, definindo incertezas $\epsilon_{i}^{}$ = x_i^o - x_i,

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 χ² de S².

Se medimos os valores I^o_i, por exemplo, o espectro, e calculamos I_i com os parâmetros x₁,..., x_k, 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}^{}$ + d²/$\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 x_i para o valor que minimizou $\tilde{{S}}$ mais um delta, digamos 5%, e minimizando novamente $\tilde{{S}}$ com este valor de x_i 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:

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 χ² 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]}$

Note que o valor médio de χ² para k graus de liberdade é k, a moda é k-2, e a variança σ²_χ²=2k.
Algumas vezes se usa um χ² 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 χ²_red seja 1, a variança continua dependente do número de graus de liberdade k

⟨χ²_red⟩=1
σ²_χ²red=2/k

Sabemos que o teorema do limite central garante que a distribuição de χ² 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 χ² é que ele assume que as incertezas das medidas são bem determinadas e não correlacionadas. Se as incertezas forem subestimadas, o valor de χ² 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 S², como por exemplo, mimizando a função discrepância, definida como:

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

e sua integral é a distribuição normal

Com esta definição

Se tivermos uma distribuição centrada em $x_o$, e de largura $\sigma^2$,

O valor da integral de x até infinito, centrada em 0 e com desvio padrão 1, é

x=0	P(>x)=0,5
x=2	P(>x)=0,02275
x=3	P(>x)=0,00135
x=4	P(>x)=0,00003167
x=5	P(>x)=2,867 × 10-7
x=6	P(>x)=9,866 × 10-10
x=7	P(>x)=1,280 × 10-12
x=8	P(>x)=6,221 × 10-16
x=9	P(>x)=1,129 × 10-19

Se a distribuição é discreta, com N valores:

de modo que precisamos normalizar nossas probabilidades

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

Distribuição gaussiana (normal) para duas variáveis não correlacionadas.

Para duas variáveis correlacionadas, x₁ medida a partir de μ₁ e x₂ medida a partir de μ₂:

onde

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 (χ²), 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