% Polynomial least-squares

function [a R2 K]=LeastSquares1(x,y,M)

N=length(x);
A=zeros(M,M);
b=zeros(M,1);

for i=1:M
    sumXY=0;
    for k=1:N
        sumXY=sumXY+y(k)*(x(k)^(i-1));
    end
    b(i)=sumXY;
    for j=1:M
        sumX=0;
        for k=1:N
            sumX=sumX+x(k)^(i+j-2);
        end
        A(i,j)=sumX;
    end
end
a=A\b;

sum=0;
for i=1:N
    sum=sum+y(i);
end
y_mean=sum/N;

sum=0;
for i=1:N
    sum=sum+(y(i)-y_mean)^2;
end
So=sum;

sum=0;
for i=1:N
    sum=sum+(y(i)-f_LS(a,x(i)))^2;
end
SE=sum;

R2=(So-SE)/So;

K=cond(A);