%% Time Varying Systems and Computation
% Klaus Diepold, Matthias Kissel, Sven Gronauer, Nehir Güzelkaya
%
%% Regression
% Regression is the problem of relating a quantity of interesty and several observed variables x_1,x_2,\dots,x_n known as the covariates, features, or independent variables.
% A regression model describes how the response is determined by the covariates through a function f applied to the features, and typically perturbed by noise z:
% y_i=f(x_i)+z_i
% We do not know this function f
% , but we have access to a dataset:

%Generation of the dataset:
n=500;                 %number of samples
sigma=0.2;             %standard deviation of the distrubition of noise
x=rand(n,1);           %x=[x_1;x_2;...;x_n]
z=sigma*randn(n,1);    %z=[z_1;z_2;...;z_n]
y=x+z;                 %x=[y_1;y_2;...;y_n]
%Plot of the dataset:
figure
scatter(x,y)
xlabel('x');
ylabel('y');
hold on

% Our ultimate goal is to predict the value y corresponding to a feature vector x that is not in the training set. To do that, we should first learn the function from the dataset, which yields an estimate \hat{f}
%  of the function f. While doing this we use the following model
% f(x_i) = a x_i+b  \quad \text{for}\: i=1,2,\dots,n.
% Now, we can convert this into a matrix notation
%  \mathbf{y}=\pmatrix{y_1  \cr \vdots \cr y_n}= \pmatrix{1 & x_1 \cr \vdots & \vdots \cr 1 & x_n}  \pmatrix{b\cr a }+  \pmatrix{z_1  \cr \vdots \cr z_n} =X \mathbf{\theta }+ \mathbf{z}

X=[ones(n,1) x];

% To find the estimation \hat{\theta} of the original coefficients {\theta}, we use least squares estimation:
% \hat{ \mathbf{\theta } } =(X^\top X)^{-1}X^\top\mathbf{y}

theta_hat=pinv(X)*y;
%Plot best fit
x_test=rand(n,1); 
y_hat=[ones(n,1) x_test]*theta_hat;
plot(x_test,y_hat,'LineWidth',3)
hold off

% Finally, we check the mean squared error of the estimation:
% ||{ \mathbf{\theta } }-\hat{ \mathbf{\theta } }||_2^2

mse=norm([0;1]-theta_hat)^2

%
% Now, let's try regression on another dataset:

%Generation of the dataset:
n=500;                
sigma=0.1;            
x=2*rand(n,1)-ones(n,1);           
z=sigma*randn(n,1);   
y=x.^2+0.5*x+z;       
%Plot of the dataset:
figure
scatter(x,y)
xlabel('x');
ylabel('y');
hold on
%Estimation
X=[ones(n,1) x x.^2];
theta_hat=pinv(X)*y;
%Plot best fit
x_test=linspace(-1,1,n)';
y_hat=[ones(n,1) x_test x_test.^2]*theta_hat;
plot(x_test,y_hat,'LineWidth',3)
hold off
% mean squared error of the estimation
mse=norm([0;0.5;1]-theta_hat)^2

%
%% Questions:
% 1. Could we use

theta_hat=inv(X'*X)*X'*y;

% instead of

theta_hat=pinv(X)*y;

% Would you expect any differences between these two commands? Explain your answer.
% 2. What would happen if number of features > n (number of training samples)? Would this change your answer to the first question?
% Hint: Think how the relation between number of features and number of training samples affects the singularity of the matrix X.