x=0:.05:1;
y= [.486 .866 .944 1.144 1.103 1.202 1.166 1.191 1.124 1.095 1.122 1.102 1.099 1.017 1.111 1.117 1.152 1.265 1.38 1.575 1.857]';

plot(x,y,'x')
title('Data for cubic Least Squares example');

%input(''); % pause for 

z= [ones(size(x)); 2*x - 1; 8*x.^2 - 8*x + 1; 32*x.^3 - 48*x.^2 + 18*x - 1]';

%coefficient matrix- symmetric
A= z'*z
b= z'*y

%find the solution
a= A\b % = (z'*z)\(z'*y)
% a =
%     1.1610
%     0.3935
%     0.0468
%     0.2396

hold on
plot (x, a(4)*(32*x.^3-48*x.^2+18*x-1) + a(3)*(8*x.^2-8*x+1) + a(2)*(2*x-1) + a(1), 'k')
% same as plot we had before: 
plot (x, 7.6687*x.^3 + -11.1282*x.^2 + 4.7259*x + 0.5747, 'r-.')

% check the error for least squares- find E_r
er= sum((y-z*a).^2)

% the coefficient of determination
r2= 1-er/sum((y-mean(y)).^2)
% therefore 97.2% of the original uncertainty has been explained by the
% nonlinear model

% but is our function well-conditioned?
cond= norm(z'*z, inf) * norm(inv(z'*z),inf)
% yes: 4.7980
log10(cond)
% therefore, expect 0 digits to be off 
% means that a small change in b will make a small change to the solution

% try it & see
b2= b+[0.01 -.01 .01 -.01]'
a2= A\b2
% a2 =
%     1.1618
%     0.3917
%     0.0483
%     0.2383
% small change in b had a small change in our solution

%Lets test using polyfit to see the best fit of the data

hold off
p1= polyfit(x,y',1);
p2=polyfit(x,y',2);
p3=polyfit(x,y',3);
plot(x,y','ko',x,polyval(p1,x), '--', x, polyval(p2, x), '-.', x, polyval(p3, x), 'LineWidth',2);
legend('data points', 'degree 1 fit', 'degree 2 fit', 'degree 3 fit');