% Skript zur Musterlösung
%
% Numerische Mathematik, SS 2013
%
% Filip Data Set Ausgleichspolynom
%
% Ingolf Busch, Juni 2013

% Daten aus filip.dat auslesen

load filip.dat

y = filip(:,1);
x = filip(:,2);
n = length(x);

% Polynomgrad festlegen
d = 10;
p = d+1;

% Bereich für den Plot
u = (-8.78:.01:-3.13)';

% NIST Referenz

beta0 = flipud([ ...
        -1467.48961422980
        -2772.17959193342
        -2316.37108160893
        -1127.97394098372
        -354.478233703349
        -75.1242017393757
        -10.8753180355343
        -1.06221498588947
        -0.670191154593408E-01
        -0.246781078275479E-02
        -0.402962525080404E-04]);
p0 = polyval(beta0,u);

f0 = polyval(beta0,x);
r0 = norm(y-f0);
rsd0 = 0.334801051324544E-02;


% Polyfit
fprintf('Polyfit\n');
beta1 = polyfit(x,y,d)';
f1 = polyval(beta1,x);
r1 = norm(y-f1);
rsd1 = r1/sqrt(n-p);
p1 = polyval(beta1,u);


% Vandermonde and backslash
fprintf('Vandermonde und backslash\n');
X = vander(x);
X = X(:,n-d:n);
beta2 = X\y;
f2 = polyval(beta2,x);
r2 = norm(y-f2);
rsd2 = r2/sqrt(n-p);
p2 = polyval(beta2,u);


% Pseudoinverse
fprintf('Pseudoinverse\n');
beta3 = pinv(X)*y;
f3 = polyval(beta3,x);
r3 = norm(y-f3);
rsd3 = r3/sqrt(n-p);
p3 = polyval(beta3,u);


% Normal equations (reject)
fprintf('Normalgleichungen(inv)\n');
beta4 = inv(X'*X)*X'*y;
f4 = polyval(beta4,x);
r4 = norm(y-f4);
rsd4 = r4/sqrt(n-p);
p4 = polyval(beta4,u);

% Normal equations (reject)
fprintf('Normalgleichungen(backslash)\n');
beta5 = (X'*X)\(X'*y);
f5 = polyval(beta5,x);
r5 = norm(y-f5);
rsd5 = r5/sqrt(n-p);
p5 = polyval(beta5,u);


% Centered
fprintf('Polyfit(centered)\n');
mu = mean(x);
sigma = std(x);
t = (x-mu)/sigma;
gamma6 = polyfit(t,y,d)';
f6 = polyval(gamma6,t);
r6 = norm(y-f5);
rsd6 =r6/sqrt(n-p);
p6 = polyval(gamma6,(u-mu)/sigma);



%Auswertungen grafisch

figure
plot(x,y,'.',u,p0)
set(legend('Data','NIST',4),'pos',[.6 .3 .18 .11])
title('NIST - Referenzlösung')
pause

figure
plot(x,y,'.',u,[p0 p1])
set(legend('Data','NIST','polyfit',4),'pos',[.6 .3 .18 .11])
title('NIST vs. polyfit')
pause

figure
plot(x,y,'.',u,[p0 p2])
set(legend('Data','NIST','backslash',4),'pos',[.6 .3 .18 .11])
title('NIST vs. backslash')
pause

figure
plot(x,y,'.',u,[p0 p3])
set(legend('Data','NIST','Pseudoinverse',4),'pos',[.6 .3 .18 .11])
title('NIST vs. Pseudoinverse')
pause

figure
plot(x,y,'.',u,[p0 p4])
set(legend('Data','NIST','Normalgleichung(inv)',4),'pos',[.6 .3 .18 .11])
title('NIST vs. Normalgleichung(inv)')
pause

figure
plot(x,y,'.',u,[p0 p5])
set(legend('Data','NIST','Normalgleichung(backslash)',4),'pos',[.6 .3 .18 .11])
title('NIST vs. Normalgleichung(backslash)')
pause

figure
plot(x,y,'.',u,[p0 p6])
set(legend('Data','NIST','polyfit(centered)',4),'pos',[.6 .3 .18 .11])
title('NIST vs. polyfit(centered)')
pause

% Auswertungen tabellarisch
fprintf('\nVergleich Residuen\n');
fprintf(1,'NIST                         %22.18f\n',r0);
fprintf(1,'polyfit                      %22.18f\n',r1);
fprintf(1,'backslash                    %22.18f\n',r2);
fprintf(1,'Pseudoinverse                %22.18f\n',r3);
fprintf(1,'Normalgleichungen(inv)       %22.18f\n',r4);
fprintf(1,'Normalgleichungen(backslash) %22.18f\n',r5);
fprintf(1,'polyfit(centered)            %22.18f\n',r6);


fprintf('\nVergleich Standardabweichung\n');
fprintf(1,'NIST                         %22.18f\n',rsd0);
fprintf(1,'polyfit                      %22.18f\n',rsd1);
fprintf(1,'backslash                    %22.18f\n',rsd2);
fprintf(1,'Pseudoinverse                %22.18f\n',rsd3);
fprintf(1,'Normalgleichungen(inv)       %22.18f\n',rsd4);
fprintf(1,'Normalgleichungen(backslash) %22.18f\n',rsd5);
fprintf(1,'polyfit(centered)            %22.18f\n',rsd6);