% make some matrices for use in exercise 2
clear
% the following gets overwritten by subsequently loading hw2data.mat,
% but you can see how I made these matrices and vectors
% A=[1  3; 5 7];
% B=[1; -2 ];
% C=[randperm(5); randperm(5);randperm(5);randperm(5);randperm(5);];
% D=reshape(primes(230),5,10);
% nt=100;
% x=[1:nt]';
% y=x./10+sin(2*pi.*x./10)+randn(nt,1);
% save hw2data.mat A B C D x y % commented out to not overwrite 

load /Users/mevans/stda/xdisk/hw2data.mat

% Problem 2
q2a=A'*B;
q2b=B'*A;
q2c=A*B;
q2d=B*B';
q2e=inv(A);
warning off MATLAB:singularMatrix % quiet error message 
q2f=inv(B*B'); % also could have asked for inv(q2d) to save cpu time
warning on MATLAB:singularMatrix  % turn warning back on

% Problem 3
q3a=C'*D;
q3b=D'*C;
q3c=C*D;
q3d=D*D';
q3e=inv(C);
q3f=inv(D*D'); % also could have asked for inv(q2d) to save cpu time
% I know these answers are correct because the operations 
% are identical to those used in Problems 2a-f; all that 
% has changed is the data upon which we operate.

% Problem 4: linear regression of y on x
% The general problem is y(t)=m1 + m2*x(t):
% y(1)=m1 + m2*x(1)
% y(2)=m1 + m2*x(2)
% y(3)=m1 + m2*x(3)
% ...
% y(nt)=m1 + m2*x(nt)
% 
% in matrix notation with suggested matrix names this is 
% Xm = y, with
nt=length(x);
X=[ones(nt,1) x];
% and solving for m gives [e.g. Sokal and Rohlf, 1995]
m=inv(X'*X)*X'*y;
% plot data and regression estimate yhat=Xm:
yhat=X*m; % best fit line through data
% calculate root-mean-squared error of the residual:
rms=sqrt((y-yhat)'*(y-yhat)./nt);[rms]
plot(x,y,'o',x,yhat,'-');grid; % plot data and regression
legend('data','linefit'); % add a legend
title('least squares regression of y on x: RMS error = 1.14 (dimensionless)'); %always title a plot
% RMS error not too big compared to sqrt(var(y))=3.16.
ylabel('y');xlabel('x') % always label axes
print -dps2 hw2fig.ps 
save hw2soln.mat