% A complete uni-variate least-squares fit
% Textbook: Part I, "Uni-Variate Linear Least-Squares", example "a complete uni-variate fit by hand"
%
% Five observations, model yhat = w1*x + w0:
%   w1 = Cov(x,y)/Var(x) = 1.7,  w0 = mean(y) - w1*mean(x) = -0.1
% The residuals sum to zero, MSE = 0.22, while the best constant model
% (yhat = mean(y)) has MSE = s_y^2 = 6.

clearvars; close all

x = [1 2 3 4 5];
y = [2 3 5 6 9];
M = length(x);

w1 = sum((x - mean(x)).*(y - mean(y)))/sum((x - mean(x)).^2);
w0 = mean(y) - w1*mean(x);
yh = w1*x + w0;
e = y - yh;

% The same fit in matrix form
X = [x' ones(M, 1)];
w = X\y';

fprintf('w1 = %.3f, w0 = %.3f (matrix form: %.3f, %.3f)\n', w1, w0, w(1), w(2))
fprintf('residuals: %s, sum = %.2g\n', mat2str(e, 3), sum(e))
fprintf('MSE = %.3f, s_y^2 = %.3f, R^2 = 1 - MSE/s_y^2 = %.4f\n', ...
    mean(e.^2), var(y, 1), 1 - mean(e.^2)/var(y, 1))
fprintf('r_xy = %.4f, r_xy^2 = %.4f\n', w1*std(x, 1)/std(y, 1), (w1*std(x, 1)/std(y, 1))^2)

figure
plot(x, y, 'ko', 'MarkerFaceColor', 'k')
hold on
plot([0.5 5.5], w1*[0.5 5.5] + w0, 'b', 'LineWidth', 1.5)
yline(mean(y), 'r--', 'LineWidth', 1.2)
for k = 1:M
    plot([x(k) x(k)], [y(k) yh(k)], 'b:', 'LineWidth', 1.2)   % residual segments
end
hold off
xlabel('x')
ylabel('y')
legend('data', sprintf('LS line, MSE = %.2f', mean(e.^2)), ...
    sprintf('constant, MSE = %.0f', var(y, 1)), 'Location', 'northwest')
grid on