% Here we fit a linear model to a noisy (simulated)
% dataset. The model is of the form
%
% y = a*x + b,
%
% And we assume given a collection of x-points and
% the corresponding y-values with some error.
%
% We will find the optimal numbers a and b in the 
% least squares sense.
%
% Samuli Siltanen November 2015

% Choose the x points
Nx   = 25;
xvec = linspace(0,1,Nx);
xvec = xvec(:);

% Pick out the a and b values used for 
% simulating data
a0 = .5;
b0 = 1;

% Compute ideal data
yvec0 = a0*xvec+b0;

% Simulate noisy data
noiselevel = 0.05;
yvec = yvec0 + noiselevel*randn(size(yvec0));

% Build the system matrix A
A = [xvec,ones(size(xvec))];

% Solve the problem in the least squares sense
tmp = inv(A.'*A)*(A.'*yvec);
a   = tmp(1);
b   = tmp(2);

% Take a look
figure(1)
clf
subplot(1,2,1)
plot(xvec,yvec0,'b')
hold on
%plot(xvec,yvec0,'k.')
plot(xvec,yvec,'r.')
axis equal

subplot(1,2,2)
plot(xvec,a*xvec+b,'r')
hold on
plot(xvec,yvec,'r.')
axis equal

