function [ a1, a0 ] = GeneralPolyRegression( x, y, J )
% Calculates coefficients for linear regression
% y = a1 x + a0

n = length(x);
a = zeros(J+1);
S = zeros(2*J+1);
Sxy = zeros(J+1);

if length(x) ~= length(y) 
    disp('Error: data vectors must be of same size');
    a = NaN;
    return
else
    fprintf('Number of data points used for linear regression: %i', n) 
end

% Create sums:
for s = 0:2*J
    S(i) = sum(x.^(i));
end

for s = 0:J
    Sxy(i) = sum(x.^(i));
end

% Create matrix:
    
Sx = sum(x);
Sy = sum(y);
Sxx = sum(x.*x);
Sxy = sum(x.*y);

a1 = (n*Sxy -Sx*Sy)/(n*Sxx - Sx.^2);
a0 = (Sxx*Sy - Sxy * Sx) / (n*Sxx - Sx.^2);

end