%linearmodel.m
%playing with linear regression (L&F pp 71-72)

%define bivariate data sample
x = [1 2 3 4 5 6]';             %column vector containing six levels of X
y = [1.5 2.7 2.7 3.5 4.5 4.8]'; %column vector containing six levels of Y

%find y-intercept (b1) and slope (b2) of linear model for this dataset
%given: y = b0 + b1*X
%hence: [b] = X\y;   using MATLAB matrix left-divide operator
X = [ones(6,1),x];   %pad X-values with 1's to represent constant slope effect
b = X\y; %left matrix division (MATLAB's shortcut for linear regression)
disp('Linear parameter estimates generated using Matrix Left Divide operation:');
disp(['   y-intercept: ',num2str(b(1))]);
disp(['         slope: ',num2str(b(2))]);
%
% Perform the same solution for the 'b' parameter vector using the
% 'long approach' without the matrix left-divide shortcut operation
% yi = b0 + b1*Xi
% y = X * b
% b = inv(X'*X) * X' * y
%   where X is a two col vector
%         first column contains all ones (represents constant intercept)
%         second column contains observed values of x (independent var)
%   y is column vector containing the y-values (dependent var)
X = [ones(6,1),x];
b = inv(X'*X)*X'*y;
disp(' ');
disp('Linear parameter estimates generated using low-level matrix algebra operations:');
disp(['   y-intercept: ',num2str(b(1))]);
disp(['         slope: ',num2str(b(2))]);

%plot best-fit linear function and data point;
xp=[0:0.01:6];
yp=b(2)*xp + b(1);
figure(1), clf;
scatter(x,y,'filled');
hold on;
axis([0 6 0 6]);
plot(xp,yp);
xlabel('X');
ylabel('Y');
str1=sprintf('slope = %04.2f', b(2));
text(1,5,str1);
str2=sprintf('y-intercept = %04.2f', b(1));
text(1,4.5,str2);

%overlay plot with y-predicted values
y2=X*b;
scatter(x,y2);
legend('Observed','Model','Predicted','Location','SouthEast')
hold off

%compute value of the discrepency function
deviations = y-y2;     %observed values - predicted values
squared_deviations = (y-y2).^2;
sum_squared_deviations = sum((y-y2).^2);
mean_squared_deviation = sum((y-y2).^2)/length(y);
rmsd = sqrt(sum((y-y2).^2)/length(y));
disp(' ');
disp(['Discrepency Function (RMSD): ',num2str(rmsd)]);
