%% Mínimos Cuadrados Totales (ejemplo Matlab)
% © 2020, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/tlsm.html , teoría en vídeo [tls1.html]
% Se tiene un muelle lineal F=k*x;

kmuelle=5;%modelo real a identificar
x=-2:.001:2; %datos limpios de entrada (inaccesibles)
N=length(x); %numero de muestras
F=kmuelle*x; %datos limpios de salida (inaccesibles)

%caso ideal análisis de componentes principales de datos limpios
datoslimpios=[F; x];
size(datoslimpios)
[Ul,Sl,Vl]=svd(datoslimpios,'econ');
Hlimpia=cov(datoslimpios')
eig(Hlimpia)
Hlimpia(1,2)/Hlimpia(2,2) %reproduce exactamente la constante del muelle.

% El análisis anterior no puede ser llevado a cabo en la práctica, porque ni x ni F "limpias" son directamente accesibles.
% Los datos para identificarlo están contaminados con ruido:

meteruido=@(v,amplitud) v+randn(size(v))*amplitud;
rng(3213)
Fmed=meteruido(F,0.35)';
xmed=meteruido(x,0.7)';

%% Mínimos Cuadrados ordinarios
% Vamos a ajustar un modelo mínimos cuadrados (LS):

th=pinv(xmed)*Fmed

% Representamos el ajuste, y el modelo real

Fest=th*xmed;
plot(xmed,Fmed,'x')

hold on
plot(x,Fmed,'.')
plot(xmed,Fest,'LineWidth',3)
plot([-2 2],[-10 10],'-.k','LineWidth',2)
hold off
legend('X Medida','X limpia', 'Modelo LS','Modelo Real')
grid on

% Los mínimos cuadrados tradicionales producen un estimado SESGADO.
%% Mínimos Cuadrados Totales
% Ahora, la versión TLS.
% Montamos matriz de datos, escalando a la misma varianza del ruido de medida.

Escalax=2;
Xi=[Fmed xmed/Escalax]';
size(Xi)
CVV=cov(Xi')
CVV(1,2)/CVV(2,2) /2 %este sería el estimado mínimos cuadrados (covarianza/varianza), el 2 es por lo del escalado.
%versión TLS con SVD (todo tiene el mismo ruido, es válido)
[U,S,V]=svd(Xi,'econ');
size(U)
size(S)
size(V)

diag(S)'/sqrt(N) %desv. típica de los componentes principales
U
-U(2,2)/U(1,2)/2 %el modelo que sale de U(1,2)*F+U(2,2)*x/2=0

% Mediante la matriz de vzas-covarianzas de la teoría:

Hmeasured=Xi*Xi'/N
[UU,D]=eig(Hmeasured);
UU
sqrt(D)

% Si restamos el mínimo valor singular por la corrupción del ruido, tendriamos los mismos vectores propios:

Hcleanaprox=Hmeasured-0.34517^2*eye(2)
[UU2,D2]=eig(Hcleanaprox);
UU2
sqrt(D2)

%
%  el modelo será

modelo=U(:,2);

% Lo escalamos para que el primer coef. sea 1:

model2=modelo/modelo(1)

% Y comparamos el ajuste:

plot(x,Fmed,'x')
hold on
plot(xmed,Fmed,'.k')
th_TLS=-model2(end)/Escalax
plot(x,[th*x' th_TLS*x'],'LineWidth',3)
hold off
grid on
legend('X limpia','Muestras','LS','TLS','Location','SouthEast')
axis tight

% El TLS no es sesgado.

