x = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50];
y = [17, 24, 31, 33, 37, 37, 40, 40, 42, 41];

% x  = [1, 2, 3, 4, 5, 6, 7];
% y = [0.5, 2.5, 2.0, 4.0, 3.5, 6.0, 5.5];

plot(x,y,'sr');
xlabel('x');
ylabel('y');

hold all;

[kcroyals, n] = size(x);

sumx = 0; sumxy = 0; st = 0;
sumy = 0; sumx2 = 0; sr = 0;

for i = 1:n
    sumx = sumx + x(i);
    sumy = sumy + y(i);
    sumxy = sumxy + x(i)*y(i);
    sumx2 = sumx2 + x(i)*x(i);
end

xm = sumx/n;
ym = sumy/n;

a1 = (n*sumxy - sumx*sumy)/(n*sumx2 - sumx*sumx);
a0 = ym - a1*xm;

for i = 1:n
    st = st + (y(i) - ym)^2;
    sr = sr + (y(i) - a1*x(i) - a0)^2;
end

syx = (sr/(n-2))^(0.5);
r2 = (st-sr)/st;

disp('Linear fit');
disp(['a0 = ' num2str(a0)]);
disp(['a1 = ' num2str(a1)]);
disp(['r^2 = ' num2str(r2)]);

JARED = a0 + a1*x;
plot(x, JARED);