%%
% Let's fit polynomials to given data

% Let's generate n data points
n=6;
% here we pick a random set of locations (x) between -1 and 1; 
%run the code many times to observe various behaviours.
x=(2*rand(n,1)-1);
x=sort(x);

y=randn(n,1);

% Now, fix the degree of the polynomial. Remember that a degree m-1
% polynomial has m coefficients a_1,..., a_m: 
% p(x)=a_1 x^(m-1)+...+a_(m-1) x + a_m

% first m=n (unique solution)
m=n;

% now, setting a=[a_1,...,a_m]^T, we have y=Va
% and since 
V=vander(x);
a=V\y;

% Let's now plot the data points (*) and the fitted polynomial over the 
% interval, say, [-1.5,1.5]

figure;
hold on
plot(x(:),y(:),'*','MarkerSize',10,'LineWidth',2)
t=min(x):0.01:max(x);
z=polyval(a,t);
plot(t,z,'LineWidth',3);

%%
% Now, let's fit with polynomials of smaller degree 

m=4;
V=vander(x);
% Pick the last m columns of V.
V=V(:,n-m+1:n);
% The command V\y will take the pseudo-inverse of V times y!
a1=V\y   % If the columns of V are not linearly independent, 
         % this gives sparsest solution. (In our case cols of V are l.i.)
a1 = pinv(V)*y   % If the columns of V are not linearly independent,
                 % this gives the solution with the smallest 2-norm.  
z1=polyval(a1,t);
plot(t,z1,'b','LineWidth',3);

%%
m=3;
V=vander(x);
V=V(:,n-m+1:n);
a2=V\y;
z2=polyval(a2,t);
plot(t,z2,'LineWidth',3);

%%
m=2;
V=vander(x);
V=V(:,n-m+1:n);
a3=V\y;
z3=polyval(a3,t);
plot(t,z3,'LineWidth',3);

%%
% Now let's experiment with structured data. Suppose the data are (noisy)
% samples of a fixed polynomial; say p(x)=x^7-x^3+1
n=8;
x=(2*rand(n,1)-1);
x=sort(x);
yp=x.^7-x.^3+1;

figure;
hold on
plot(x(:),yp(:),'b*','MarkerSize',10,'LineWidth',2)


t=min(x):0.01:max(x);



%%
% Suppose now you don't know p(x) and you want fit the data using the
% approach above. Use m=8 (as that is the right degree). Let's see what
% happens. 

m=8;
V=vander(x);
V=V(:,n-m+1:n);
a=V\yp;
z=polyval(a,t);
plot(t,z,'b','LineWidth',3);

%%

% Now, add some small noise and fit again

ypn=yp+0.1*rand(size(yp));

m=8;
V=vander(x);
V=V(:,n-m+1:n);
a=V\ypn;
z1=polyval(a,t);
plot(x(:),ypn(:),'r*','MarkerSize',10,'LineWidth',2)

plot(t,z1,'r','LineWidth',3);

%%
% Let's fit with a polynomial of smaller degree

m=6

figure;
hold on
plot(x(:),yp(:),'b*','MarkerSize',10,'LineWidth',2)

V=vander(x);
V=V(:,n-m+1:n);
a=V\yp;
z=polyval(a,t);
plot(t,z,'b','LineWidth',3);

a=V\ypn;
z1=polyval(a,t);
plot(x(:),ypn(:),'r*','MarkerSize',10,'LineWidth',2)

plot(t,z1,'r','LineWidth',3);

%%
% It all depends on how many data points we have. Say we have way more...


n=50;
x=(2*rand(n,1)-1);
x=sort(x);
yp=x.^7-x.^3+1;

t=min(x):0.01:max(x);

figure;
hold on
plot(x(:),yp(:),'b*','MarkerSize',10,'LineWidth',2)

% Suppose now you don't know p(x) and you want fit the data using the
% approach above. Use m=8 (as that is the right degree). Let's see what
% happens. 

m=8;
V=vander(x);
V=V(:,n-m+1:n);
a=V\yp;
z=polyval(a,t);
plot(t,z,'b','LineWidth',3);

%%
% Now, add some small noise and fit again

ypn=yp+0.1*rand(size(yp));

m=8;
V=vander(x);
V=V(:,n-m+1:n);
a=V\ypn;
z1=polyval(a,t);
plot(x(:),ypn(:),'r*','MarkerSize',10,'LineWidth',2)

plot(t,z1,'r','LineWidth',3);





