O método dos mínimos quadrados é uma das ideias mais frequentemente usadas na ciência. A primeira referência ao método aparece como apêndice do livro “Nouvelles méthodes pour détermination des orbites de comètes” publicado por Legendre em 1805. Mesmo assim na época Gauss revindicou a prioridade da descoberta, alegando que vinha usando o método desde 1795. Plackett faz um resgate detalhado desta história.
Considere que temos um conjunto de n medidas experimentais da quantidade y para diferentes valores de uma variável experimental x, por exemplo, a absorbância em função da concentração. Considere também existe um modelo teórico que propõe uma relação entre y e x, tal que y=f(x;{pj}j) onde {pj}j são um conjunto de m parâmetros do modelo. Por exemplo, a lei de Beer estabelece que a absorbância da espécie X no comprimento de onda l é dado por Al =ell[X]. Assim se y=Al e x=[X] então y=f(x;{pj}j)=p1+p2x, onde p1=0 e p2=ell.
Normalmente os parâmetros {pi}i do modelo f(x;{pi}i) são desconhecidos e o nosso objetivo é determina-los a partir do dados experimentais {xi,yi}i, com i=1…n.
A função y=p1+p2x fornece no plano xy os pontos geométricos da reta, visto que Dy=p2Dx então cada incremento em x de Dx leva ao mesmo incremento Dy em y. Outra forma de ver isto é que dy/dx para qualquer x tem sempre o mesmo valor p2. O parâmetro p1 corresponde ao valor de y quando x=0, isto é y=p1 quando x=0. O parâmetro p2 corresponde a inclinação da reta (que é a mesma para todo x), já que dy/dx=p2=tan(q) onde q é o ângulo entre a reta e o eixo x. Portanto, para cada valor de p1 e p2 temos uma reta diferente. A questão então é: Qual o valor de p1 e p2 da reta que mais se aproxima dos dados experimentais {xi,yi}i=1n.
O problema acima como colocado está vago e impreciso: o que significa “que mais se aproxima dos dados experimentais”? Vamos considerar que foram calculados o valor previsto da quantidade y, usando o modelo yical=f(x;{pi}i)=p1+p2xi para os n valores experimentais da variável x. A diferença entre o valor experimental e calculado para cada x é chamado de ei, assim
ei=yical-yi
Agora podemos definir o que chamamos de função erro, como sendo o somatório dos erros ao quadrado, isto é
E({pj}j)=Si ei2
Veja que E é uma função dos parâmetros do modelo, isto é para cada conjunto {pj}j temos diferentes valores para yical, portanto, diferentes valores para ei e, portanto, um valor diferente para E. Este é o critério que vamos usar então para definir a função que mais se aproxima dos dados experimentais. Neste caso, a melhor função (ou conjunto de parâmetro {pj}j para uma função expecífica) é aquela que minimiza a função E*. Podemos escrever isto como
{pj}jótimo = min_{{pj}j} E({pj}j)
*Dizemos função E quando depende dos parâmetros {pj}j para um f específico e dizemos funcional E quando depende de f. Por hora vamos considerar que E depende de um conjunto de parâmetros.
Vamos considerar o problema numa perspectiva diferente agora. Considere o ponto P(xi,yical), tal que yical=f(x;{pi}i). agora considere um ponto externo xp
% regressão ortogonal, jan 2023, nhtl, unifal
function []=ortoregr()
clear all
close all
global x y
global N
rng(899823993)
x1=1:10; x=x1+rand(size(x1));
y1=2*x1+3; y=y1+rand(size(x1));
N=length(x);
p0=[1.2, 2.1];
options = optimset('PlotFcns',@optimplotfval,'TolFun',1e-6,'TolX',1e-6);
p = fminsearch(@Custo,p0,options)
pc=linortfit(x,y); [pc(2) pc(1)]
pm = polyfit(x,y,1)
xgraf=x(1)-1:.1:x(end)+1;
figure(2)
plot(x,y,'*r')
hold on
plot(xgraf,2*xgraf+3,'r-')
plot(xgraf,p(1)*xgraf+p(2),'k-')
plot(xgraf,pc(2)*xgraf+pc(1),'b-')
plot(xgraf,pm(1)*xgraf+pm(2),'g-')
end
function C=Custo(p)
global x y
global N
C=0;
for i=1:N
xp=minE(i,p);
d2=(y(i)-func(xp,p))^2+(x(i)-xp)^2;
C=C+d2;
end
end
function xp=minE(i,p)
global x y
xp = fzero(@(xp) Energia(xp,i,p),x(i));
end
function E=Energia(xp,i,p)
global x y
E=x(i)-xp+dfunc(xp,p)*(y(i)-func(xp,p));
end
function f=func(xp,p)
a=p(1); b=p(2);
f=a*xp+b;
end
function df=dfunc(xp,p)
a=p(1);
df=a;
end
