Algoritmo en MATLAB para la aproximación lineal por el método de los mínimos cuadrados.
Por: Carlos Armando De Castro P.
El siguiente algoritmo recibe un número arbitrario de pares de datos en la forma de una matriz de 2*n, donde las abcisas se encuentran en la primera fila (o renglón) y las ordenadas en la segunda fila de la matriz, devolviendo la pendiente m y el intercepto b de la recta que interpola a los datos y además entrega su gráfica:
function [m,b]=mincuadlin(X) n=length(X(1,:)); A=0; B=0; C=0; D=0;
for i=1:n; A=A+X(1,i); B=B+X(2,i); C=C+(X(1,i))^2; D=D+X(1,i)*X(2,i); end
m=(n*D-A*B)/(n*C-A^2); b=(C*B-D*A)/(n*C-A^2);
for i=1:n; hold on; plot (X(1,i),X(2,i),'*','MarkerEd (X(1,i),X(2,i),'*','MarkerEdgeColor','r','LineWidth',1); geColor','r','LineWidth',1); end
x=X(1,1):1:X(1,n); y=m*x+b; plot(x,y,'b'); title('Aproximación lineal por mínimos cuadrados.'); cuadrados.') ;
Por ejemplo, para los datos {(1,0),(2,3),(3,4),(4,-6),(5,2),(6,4),(7,0),(8,4) {(1,0),(2,3),(3,4),(4,-6),(5,2),(6,4),(7,0),(8,4),(9,3)}, ,(9,3)}, se escribe en el Command Window :
>>X=[1 2 3 4 5 6 7 8 9; 0 3 4 -6 2 4 0 4 3]; >>[m,b]=mincuadlin(X) Y el programa entrega los resultados: m = 0.2833 b = 0.1389
Descripción: Este código realiza ajustes varios de regresión por mínimos cuadrados. Los ajustes son exponencial, de potencia, logarítmico y polinomial. Al final, se visualiza en pantalla la ecuación del modelo ajustado y su gráfica correspondiente.
% Regresion minimos cuadrados % Ajustes varios % Realizado por Ing E.Porto % El programa realiza varios ajustes: exponencial, de potencias, % logarítmico y polinomial
clc;clear; %Ingreso de datos. Es lo único que tienes que digitar. x=[1 2 3 4];y=[0 0.5 2 4]; %Digita aqui los datos de las parejas X y Y
%Opciones de inicio del programa n=length(y); opcion=0;sg=1;ajuste=1;xn=x;yn=y;
%Menú para escoger análisis de regresión while opcion==0 k=menu('Escoja analisis regresion','A*exp(B*X)','A*B^X','A*X^B','A+B*Ln(X)','Polinomio'); if k==1 opcion=1; m=1; if any(y<0) ajuste=0; else yn=log(y); end elseif k==2 opcion=1; m=1;
if any(y<0) ajuste=0; else yn=log(y); end elseif k==3 opcion=1; m=1; if any(y<0)|any(x<0) ajuste=0; else yn=log(y); xn=log(x); end elseif k==4 opcion=1; m=1; if any(x<0) ajuste=0; else xn=log(x); end elseif k==5 opcion=1; m=input(' Grado del polinomio = '); if m>n ajuste=0; end end end
% Cuerpo del programa if ajuste==0 fprintf(' No se puede realizar el ajuste ');
else for i=1:m+1 for j=1:m+1 sx(i,j)=sum(xn.^(i+j-2)); end sy(i)=sum(yn.*xn.^(i-1)); end % Presentacion de resultados fprintf(' Matriz de sumatorias '); disp([sx sy']); c=sxsy'; xx=linspace(min(x),max(x)); fprintf(' Curva ajustada: '); if k==1 fprintf(' Y = %g * exp(%g * X) n',exp(c(1)),c(2)); yy=exp(c(1))*exp(c(2)*xx); ya=exp(c(1))*exp(c(2)*x); elseif k==2 fprintf(' Y = %g * %g ^ X n',exp(c(1)),exp(c(2))); yy=exp(c(1))*exp(c(2)).^xx; ya=exp(c(1))*exp(c(2)).^x; elseif k==3 fprintf(' Y = %g * X ^ %g n',exp(c(1)),c(2)); yy=exp(c(1))*xx.^c(2); ya=exp(c(1))*x.^c(2); elseif k==4 fprintf(' Y = %g + %g * LnX n',c(1),c(2)); yy=c(1)+c(2)*log(xx); ya=c(1)+c(2)*log(x); elseif k==5 for w=1:m+1 if c(w)<0 sg='-'; else
sg='+'; end fprintf('%s %g X^%g ',sg,abs(c(w)),w-1); end cn=flipud(c); ya=polyval(cn,x); yy=polyval(cn,xx); end
%Impresión en pantalla de resultados fprintf(' Presione cualquier tecla para ver la grafica del ajuste '); pause %Gráfica del modelo ajustado y de los puntos originales plot(xx,yy,x,y,'.r'); pause
end
