% This is a small script to help you review some of the
% concepts covered during class for the linear regression
% lecture

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% First let's set up some data. In this example the data comes from
% a simple line with slope .75 and Y-axis intercept=2

w=.75;
b=2;

% Let's get a few locations randomly along this line between x=0 and x=5;

x=5*rand(15,1);   % Get 15 random x locations between 0 and 5
	
t=(w*x)+b;	  % Compute the exact value of t from the line equation

% Plot the noiseless data

figure(1);clf;
plot(x,t,'bo');hold on;
title('Data sampled from a line. Blue is noiseless, magenta is with noise');

% Looks nice, but unfortunately we don't generally work with perfect
% data, sources of error: Measurement problems, transmission distortion,
% random noise or influences from the environment.

%% We will look at typical noise models later on, in particular, we will
% use Gaussian functions to model nose quite often. Let's for now
% use Matlab to generate some random Gaussian noise, with mean of
% zero and a standard deviation of .5 

t_noisy=t+(.5*randn(length(t),1));

% Plot it on the same figure in red to see what it looks like

plot(x,t_noisy,'ms');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, let's have some fun with Matlab:
%
% Go to the Tools menu, and select Basic Fitting
% This will open up a menu
% - Select the data source corresponding to the BLUE points first
%    (ie data set 1)
% - Then check on the 'linear' checkbox
%   You will see the best-fit line estimated from our samples
%   using the same procedure we covered in lecture
% - Try fitting the data with other functions
%   
% Question: Do all the functions approximate the data equally well?
%           what is going on here?
%           Is there a point at which Matlab complains about something?
%          
% - Now, select the data source corresponding to the MAGENTA points
% - Try the linear fit again
% - With the linear fit on, try other functions
%
% Question: Are the results similar to the noiseless case? can all
%           functions approximate the data equally well?
%           if not, which function does a better job?
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Clear this plot and let's get some more data and more noise

clf;

% Now let's see what happens when we have more data, and also, more noise!
% Let's take twice the amount of data

x=5*randn(30,1);
t=(w*x)+b;
t_noisy_1=t+(.1*randn(length(t),1));	% Random noise stDev=.1
t_noisy_2=t+(.75*randn(length(t),1));	% Random noise stDev=.75

figure(1);plot(x,t_noisy_1,'bs');hold on;
plot(x,t_noisy_2,'m*');
title('Noisy data with different amounts of noise');

% Try out fitting these data with different functions
% like above... do you get similar results? overall,
% what would you say is the function that best 
% approximates the data?


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Least Squares Regression 
%
% We now know Matlab can solve this problem, but let's do it
% ourselves using the Least Squares closed-form solution for
% parameters w and b that we discussed 
%
% Remember, the exact paramaters for our line are
% w=.75, b=2
%
% From the notes: w=[sum_i (t_i-T)(x_i-X)]/sum_i (x_i-X)^2
%                 b=T-wX
%
% where T is the average t value for the input samples
% and X is the average x value for the input samples
%

% Let's look at the noiseless case first (t and x in our workspace)

X_noiseless=mean(x);	% Compute average x
T_noiseless=mean(t);	% Compute average t

w_noiseless=sum((t-T_noiseless).*(x-X_noiseless))/sum((x-X_noiseless).^2)
b_noiseless=T_noiseless-(w_noiseless*X_noiseless)

% We get the exact solution, which is as we expected since in the noiseless
% case the optimal solution indeed corresponds to the original model.
    
% Now see what happens with the noisy cases, first, with little noise
% (t_noisy1 in our workspace)

% The mean x value doesn't change! notice that we only added noise to
% the t measurements. Hence, we don't need to recompute the averade
% x value and we simply use X_noiseless.

T_noisy1=mean(t_noisy_1);	% Mean t for noisy 1 data

w_noisy1=sum((t_noisy_1-T_noisy1).*(x-X_noiseless))/sum((x-X_noiseless).^2)
b_noisy1=T_noisy1-(w_noisy1*X_noiseless)

% Close! but the effect of noise is clear, see what happens with
% t_noisy_2

T_noisy2=mean(t_noisy_2);	% Mean t for noisy 2 data

w_noisy2=sum((t_noisy_2-T_noisy2).*(x-X_noiseless))/sum((x-X_noiseless).^2)
b_noisy2=T_noisy2-(w_noisy2*X_noiseless)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now let's look at multivariate regression
% 
% Let's simulate some data from the model t = w1 + w2*x2 + w3*x3 + w4*x4
% with the parameters w = [w1 w2 w3 w4]

clear;
w = [1 .4 .2 -.6]';

% Create 50 random data points with x2 %in% [0;5], x3 %in% [0;10], 
% x4 %in% [0;8]
N = 50;
X = [ones(N, 1) rand(N, 3)];
X = repmat([1, 5, 10, 8], size(X, 1), 1) .* X;

s = .3; % Standard deviation
t = X * w + s * randn(N, 1);

% First let's do it the WRONG way!
% Solving the normal equations directly
%
% In this little toy example, inverting the matrix directly does no harm,
% but in general, this is very unstable numerically.
w_hat1 = inv(X'*X)*X'*t;

% Instead, use the magic 'backslash' operator
% Easy, fast and numerically stable (using the pseudo-inverse of X)
w_hat2 = X \ t; 

% Compute the fitted values
t_hat = X*w_hat2;

% What's the error?
E = norm(t - t_hat);
