Approximation by polynomials using the least squares method

Using Matlab function polyfit

>> X =  [-1 1 2 -1 0 2];   % x-coordinates of the points
>> Y =  [ 8 4 5 7 4 6];    % y-coordinates of the points

>> deg = 2;                % degree of the polynomial
>> p = polyfit(X, Y, deg)  % computation of coeff. of the polynomial

>> step = 0.2;             % step-size of the graph of the polynomial
>> xx = min(X):step:max(X);  % division on x-axis 
>> yy = polyval(p, xx);    % corresponding values of the polynomial
>> plot(X, Y, 'ro', xx, yy);

>> Yp = polyval(p, X);     % approximate values at given points
>> D = abs(Y-Yp);          % distances from given values
>> delta = sqrt(sum(D.^2)) % standard deviation
>> max_dev = max(D)        % maximal deviation

Experiment with changing the degree of the polynomial from constant (i.e. average value, degree 0), a line (degree 1), up to interpolation polynomial (it can be obtained for distinct x-coordinates only).

Using a projection

You can use basics of Matlab instead of the function polyfit:

>> X =  [-1 1 2 -1 0 2];  % x-coordinates of the points
>> Y =  [ 8 4 5 7 4 6];   % y-coordinates of the points

>> deg = 2;               % degree of the polynomial

>>   %%% computation of coeff. of the polynomial:
>> n = length(X); 
>> Q = ones(n,1);        % the first column - ones
>> for k = 1:deg
>>   Q = [Q , (X.^k)' ]; % for every power k, add one column
>> end

>>   %%%   Q' performs a projection to lin. space
>>   %%%   generated by columns of Q:
>> A = Q' * Q;          
>> b = Q' * Y'; 

>> q = A\b;     % computation of coeff. of the polynomial

>>   %%%  graph of the polynomial
>> p = flip(q)  % reversing the elements in q        
>> step = 0.2;          % step-size of the graph of the polynomial
>> min(X):step:max(X);  % division on x-axis 
>> yy = polyval(p, xx);  % values of the polynomial at xx
>> plot(X, Y, 'ro', xx, yy);

>>   %%% checking the quality of approximation
>> Yp = polyval(p, X);  % approximate values at given points
>> D = abs(Y-Yp);       % distances from given values
>> delta = sqrt( sum( D.^2) ) % standard deviation
>> max_dev = max(D)     % maximal deviation