% Listing 2.1 of the book
% For Example 2.3, Page 34

% Arun K. Tangirala
% March 21, 2015
% Freely distributed under no warranty

%% Create input-output data
uk = rand(200,1); yk = 5*uk + 2;

% Add noise to get the measurement
ek = randn(200,1);
ykm1 = yk + ek*sqrt(var(yk)/100);	% SNR 100
ykm2 = yk + ek*sqrt(var(yk)/10);	% SNR 10

%% Plot data and use curve fitting to obtain the best fits
figure; plot(uk,ykm1,'x','linewidth',1);
figure; plot(uk,ykm2,'x','linewidth',1);

%% Compute LS estimates and their standard errors
Phi = [uk ones(200,1)]
thetah1 = Phi \ ykm1;
err1 = ykm1 - Phi*thetah1;
errth1 = sqrt(diag(inv(Phi'*Phi)*sum(err1.^2)/198));

thetah2 = Phi \ ykm2;
err2 = ykm2 - Phi*thetah2;
errth2 = sqrt(diag(inv(Phi'*Phi)*sum(err2.^2)/198));