% 
% confronto tra algoritmi per i minimi quadrati discreti:
%
coord = (1:50);
gradopol = 5;  % "gradopol": grado del polinomio che genera i dati osservati:
% "x_veri": valori veri dei coefficienti del polinomio che genera i dati osservati:
x_veri = [2.987654321 5.987654321 4.987654321 1.987654321 3.987654321 0.987654321 2.987654321 1.987654321]';
x_veri = x_veri(end-(gradopol:-1:0));
% "A": matrice del modello lineare
A = vander(coord);   A = A(:,end-(gradopol:-1:0));
%
% caso di misure perfette :
b = polyval(x_veri,coord)';
[Qh,Rh]=qr(A);
x1 = Rh \ (Qh'*b);
[Qm,Rm]=qr_MGS(A);
x2 = Rm \ (Qm'*b);
x3 = (A'*A) \ (A'*b);
disp('risultato con misure perfette:')
[(x1-x_veri)./x_veri  (x2-x_veri)./x_veri  (x3-x_veri)./x_veri]
disp('QR_Householder è nettamente il migliore !')
pause
disp(' '); disp('-----------------------------------------------'); disp(' ');
% caso di misure rumorose :
stderr = 1e-3;
br = polyval(x_veri,coord)' + stderr*randn(length(coord),1);
[Qhr,Rhr]=qr(A);
x1r = Rhr \ (Qhr'*br);
[Qmr,Rmr]=qr_MGS(A);
x2r = Rmr \ (Qmr'*br);
x3r = (A'*A) \ (A'*br);
disp(['risultato con misure rumorose: (stderr = ' num2str(stderr) ')'])
[(x1r-x_veri)./x_veri  (x2r-x_veri)./x_veri  (x3r-x_veri)./x_veri]
disp('il rumore nelle misure corrompe la precisione degli algoritmi:')
disp('per effetto del rumore nelle misure, i tre algoritmi sono adesso equivalenti !')
pause
disp(' ');
disp(' ');
disp(' '); disp('-----------------------------------------------'); disp(' ');

%
% curiosità: vediamo che risultati darebbe l'interpolazione :
%
% su misure perfette (primi campioni):
disp('interpolazione con misure perfette (primi campioni) :')
(A(1:gradopol+1,:)\b(1:gradopol+1) - x_veri)./x_veri
pause
disp(' '); disp('-----------------------------------------------'); disp(' ');
%
% su misure perfette (ultimi campioni):
disp('interpolazione con misure perfette (ultimi campioni) :')
(A(end-(gradopol:-1:0),:)\b(end-(gradopol:-1:0)) - x_veri)./x_veri
pause
% vediamo perchè questa differenza :
disp(['cond(A(1:gradopol+1,:)) = ' num2str(cond(A(1:gradopol+1,:)))]);
disp(['cond(A(end-(gradopol:-1:0),:)) = ' num2str(cond(A(end-(gradopol:-1:0),:)))]);
pause
disp(' ');
disp(' ');
disp(' '); disp('-----------------------------------------------'); disp(' ');
%
% su misure rumorose :
disp(['interpolazione (primi campioni) con misure rumorose: (stderr = ' num2str(stderr) ')'])
(A(1:gradopol+1,:)\br(1:gradopol+1) - x_veri)./x_veri
disp('l''interpolazione è più sensibile al rumore di quanto lo siano i minimi quadrati!')
