#Métodos Numéricos - 2023. IMERL, Fing
#Código de ejemplo para el ejercicio 3 del práctico 5
t = [1:25];
y = [5.0291, 6.5099, 5.3666, 4.1272, 4.2948, 6.1261, 12.5140, 10.0502, 9.1614, 7.5677, 7.2920, 10.0357, 11.0708, 13.4045,12.8415, 11.9666, 11.0765, 11.7774, 14.5701, 17.0440, 17.0398, 15.9069, 15.4850, 15.5112, 17.6572];

m=25;

#parteB=1 si se quiere ejecutar el código correspondiente.
#Idem para parteC.
parteB=0;
parteC=1;

##################################################
#Parte a) Ajuste lineal.
##################################################
#Calcular coeficientes del ajuste lineal
A=[ones(m,1),t'];
x=(A'*A)\(A'*y');

#Calcular valores del modelo en el vector t
yModelo=x(1)+x(2)*t;

#Graficar puntos y recta
figure(1);
scatter(t,y,50,'r','s');
hold on;
grid on;
set(gca,'fontsize', 26);
title('Datos vs. Ajuste lineal','FontSize',28);
plot([1:0.01:25],x(1)+x(2)*[1:0.01:25],'LineWidth',2,'k');
legend('Datos','Ajuste lineal',"location", 'northwest');

#Graficar residuo
figure(2)
scatter(t,abs(y-yModelo),50,'b','s');
hold on;
plot(t,abs(y-yModelo),'LineWidth',2,'b');
title('Residuo','FontSize',28);
grid on;
set(gca,'fontsize', 26);

#Hallo el índice donde se da el outlier
[resOut,indOut]=max(abs(y-yModelo));
indOut

#Conclusión: el dato con t=indOut es un outlier.
t_=t;
t_(indOut)=[];
y_=y;
y_(indOut)=[];

##################################################
#Parte b) Eliminar el outlier y rehacer el ajuste.
##################################################
if parteB==1
#Calcular coeficientes del ajuste lineal
A=[ones(m-1,1),t_'];
x=(A'*A)\(A'*y_');

#Calcular valores del modelo en el vector t
yModelo_=x(1)+x(2)*t_;

#Graficar puntos y recta
figure(3);
scatter(t_,y_,50,'r','s');
hold on;
grid on;
set(gca,'fontsize', 26);
title('Datos vs. Ajuste lineal - Sin outlier','FontSize',28);
plot([1:0.01:25],x(1)+x(2)*[1:0.01:25],'LineWidth',2,'k');
legend('Datos','Ajuste lineal',"location", 'northwest');

#Graficar residuo
figure(4)
plot(t_,abs(y_-yModelo_),'LineWidth',2,'b');
hold on;
scatter(t_,abs(y_-yModelo_),50,'b','s');
title('Residuo - Sin outlier','FontSize',28);
grid on;
set(gca,'fontsize', 26);

endif
##################################################
#Parte c) Ajustar otro modelo, teniendo en cuenta
#el patrón de los residuos de la parte b)
##################################################
#Calcular coeficientes del modelo
if parteC==1

AC=[ones(m-1,1),t_',sin(t_)'];
xC=(AC'*AC)\(AC'*y_');

#Calcular valores del modelo en el vector t
yModeloC=xC(1)+xC(2)*t_+xC(3)*sin(t_);

#Graficar puntos y modelo
figure(5);
scatter(t_,y_,50,'r','s');
hold on;
grid on;
set(gca,'fontsize', 26);
title('Datos vs. Ajuste modelo C - Sin outlier','FontSize',28);
plot([1:0.01:25],xC(1)+xC(2)*[1:0.01:25]+xC(3)*sin([1:0.01:25]),'LineWidth',2,'k');
legend('Datos','Ajuste lineal',"location", 'northwest');
endif