%general algorithm for least squares approximation

m = 100;               %m = number of data points
n = 4;                 %n-1 = degree of polynomial approximation

a = linspace(-2,2,m);  %(a(i),b(i)) are the data points
b = zeros(m,1);        %here they are generated artificially

for i=1:m
    b(i) = a(i)^3 + 4*(rand()-0.5);    %rand() = random number in (0,1)
end

%in this example, b(i) = a(i)^3 + random number in (-2,2)
%so the polynomial approximation should be close to P(x) = x^3

for i=1:m
    for j=1:n
        A(i,j) = a(i)^(j-1);
    end
end

A_t = transpose(A);
c = linsolve(A_t*A,A_t*b)          %c = coefficients of best fit polynomial

plot(a,b,'or')                     %plots data points

x = linspace(min(a),max(a),1000);
y = zeros(1,1000);                 %y = best fit function evaluated on x

for i=1:1000
    y(i) = c(1);
    r = 1;
    for j=2:n
        r = r*x(i);
        y(i) = y(i) + c(j)*r;
    end
end

plot(a,b,'or',x,y,'-b')        %plots points and polynomial approximation
