% Ajuste por minimos cuadrados en una base generica
clear all
close all
x =[0.16000.97000.96000.49000.80000.14000.42000.92000.79000.9600]
y =[0.52080.75830.73960.34140.50110.53680.35490.66930.48990.7396]
f1=inline('ones(size(x))');
f2=inline('x.^2');
f3=inline('sin(x)');
f4=inline('log10(x)');
X=[f1(x);f2(x);f3(x);f4(x)]';
Y=y';
s=(X'*X)\(X'*Y)
norm(c(1)*f1(x)+c(2)*f2(x)+c(3)*f3(x)+c(4)*f4(x)-y)

xx=linspace(01);
yy=c(1)*f1(xx)+c(2)*f2(xx)+c(3)*f3(xx)+c(4)*f4(xx);
plot(xy'*'xxyy);grid on;
xlabel('abcisas');ylabel('ordenadas');
legend('Datos''Ajuste')
title('Minimos cuadrados')
%% Ejemplos
%
% Nº 1
clear all
close all
x=[1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, ...
   1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008];
y=[1112186, 1115016, 1123301, 1125419, 1128372, 1093937, 1099296, 1112415,...
   1117762, 1117370, 1087885, 1081834, 1084314, 1076567, 1075329, 1073971,...
   1075381, 1073761, 1076635, 1076896, 1074862, 1080138];
f1=inline('ones(size(x))');
f2=inline('1./x');
X=[f1(x);f2(x)]';
Y=y';
s1=(X'*X)\(X'*Y)
e1=norm(s1(1)*f1(x)+s1(2)*f2(x)-y)
xx=linspace(1986,2008);
yy1=s1(1)*f1(xx)+s1(2)*f2(xx);
%
f3=inline('ones(size(x))');
f4=inline('x');
f5=inline('sin(x)');
X=[f3(x);f4(x);f5(x)]';
Y=y';
s2=(X'*X)\(X'*Y)
e2=norm(s2(1)*f1(x)+s2(2)*f2(x)+s2(3)*f3(x)-y)
yy2=s2(1)*f3(xx)+s2(2)*f4(xx)+s2(3)*f5(xx);
%
plot(x,y,'*',xx,yy1,xx,yy2);grid on;
xlabel('año');ylabel('poblacion');
legend('Datos','Ajuste1','Ajuste2')
title('Poblacion en Asturias')
%
% Nº 2
vary=var(y);N=length(x);
for i=1:4;
    [p,S,mu] = polyfit(x,y,i)[yy,delta] = polyval(p,xx,S,mu);
    subplot(2,2,i);plot(x,y,'*',xx,yy,'r',xx,yy-delta,'r:',xx,yy+delta,'r:');
    grid on;xlabel('año');ylabel('poblacion');
    title(['Poblacion en Asturias (grado ',num2str(i),')']);
    text(0.8,0.9,['error: ',num2str(S.normr)],'Units','Normalized');
    r2=1-S.normr.^2/vary/(N-i);
    text(0.8,0.8,['r: ',num2str(sqrt(r2))],'Units','Normalized');
end
%
% Nº 3 Representación de funciones
f=inline('a*x.^b');txt={};
a=[1,2,3];la=length(a);
b=[-1 2];lb=length(b);
xx=linspace(1,3,51);
yy=zeros(la*lb,length(xx));
for ia=1:la;
    for ib=1:lb;
        yy(ib+(ia-1)*lb,:)=f(a(ia),b(ib),xx);
        txt=[txt,['[a,b]=[',num2str(a(ia)),',',num2str(b(ib)),']']];
    end
end
%semilogy(xx,yy);grid on
plot(xx,yy);grid on
legend(txt,'Location','Best')
title(char(f))
%
f=inline('a*exp(x.*b)');txt={};
a=[1,2,3];la=length(a);
b=[-1 2];lb=length(b);
xx=linspace(0,1,51);
yy=zeros(la*lb,length(xx));
for ia=1:la;
    for ib=1:lb;
        yy(ib+(ia-1)*lb,:)=f(a(ia),b(ib),xx);
        txt=[txt,['[a,b]=[',num2str(a(ia)),',',num2str(b(ib)),']']];
    end
end
%semilogy(xx,yy);grid on
plot(xx,yy);grid on
legend(txt,'Location','Best')
title(char(f))
%
% Ajuste no lineal - linealizable
%
x=1:5;f=inline('log(x+2).^2');y=f(x);
% Potencial
p1 = polyfit(log(x),log(y),1) 
b=p1(1),a=exp(p1(2))
ya = a*x.^b; % dos formas de evaluar la función de ajuste
ya1 = exp(polyval(p1,log(x)));
er1 = norm(y-ya)
% Exponencial
p2 = polyfit(x,log(y),1)
b=p2(1),a=exp(p2(2))
ya = a*exp(x*b); % dos formas de evaluar la función de ajuste
ya2 = exp(polyval(p2,x));
er2 = norm(y-ya)
% Racional 
p3 = polyfit(1./x,1./y,1)
ya = x./(p3(2)*x+p3(1)); % dos formas de evaluar la función de ajuste
ya3 = 1./polyval(p3,1./x);
er3 = norm(y-ya)
% Comparación
plot(x,[y;ya1;ya2;ya3],'o');hold on
xx=linspace(x(1),x(end));
plot(xx,[f(xx);exp(polyval(p1,log(xx)));exp(polyval(p2,xx));1./polyval(p3,1./xx)]);
legend('Datos','Potecial','Exponencial','Racional','Location','Best');hold off
grid on


%
% Representación de las funciones para diferentes valores de los parámetros
f=inline('1./(a*x+b)');txt={};
a=[1,2,3];la=length(a);
b=[0 2];lb=length(b);
xx=linspace(0.5,2,51);
yy=zeros(la*lb,length(xx));
for ia=1:la;
    for ib=1:lb;
        yy(ib+(ia-1)*lb,:)=f(a(ia),b(ib),xx);
        txt=[txt,['[a,b]=[',num2str(a(ia)),',',num2str(b(ib)),']']];
    end
end
%semilogy(xx,yy);grid on
plot(xx,yy);grid on
legend(txt,'Location','Best')
title(char(f))
%
% Nº 4 Lectura del fichero de datos
% Se genera 'data', 'textdata' y 'rowhead'
uiimport('IG_pr10D1.txt')
x=data(:,1);y=data(:,2);
% No vale log(x) al ser x=0
p = polyfit(x,log(y),1)
e = norm(y-exp(polyval(p2,x)))