clc
close all
clear

directory = '/home/victor/Dropbox/Latex_Reports/Proyecto_Integrador_Señales_y_Sistemas/Linear_Regression/';
filename = 'Advertising.csv';

M = csvread(filename,1,1);

% TV	radio	newspaper	sales
TV = M(:,1);
radio = M(:,2);
newspaper = M(:,3);
sales = M(:,4);

figure(1)
subplot 311
plot(TV,sales,'o','markersize',4)
xlabel('TV')
ylabel('sales')

subplot 312
plot(radio,sales,'o','markersize',4)
xlabel('radio')
ylabel('sales')

subplot 313
plot(newspaper,sales,'o','markersize',4)
xlabel('newspapers')
ylabel('sales')
%% one predictor

% TV
X = [ones(length(sales),1) TV];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

subplot 311
hold on
plot(TV,salesPredict,'r','linewidth',2)
hold off

R2 = stats(1);
Fstat = stats(2);
pValue = stats(3);
var = stats(4);

text(10,25,['R^2 = ' num2str(R2)])
text(10,19,['y = ' num2str(b(2),3) ' TV + ' num2str(b(1),3)])

% radio
X = [ones(length(sales),1) radio];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

subplot 312
hold on
plot(radio,salesPredict,'r','linewidth',2)
hold off

R2 = stats(1);
Fstat = stats(2);
pValue = stats(3);
var = stats(4);

text(2,25,['R^2 = ' num2str(R2)])
text(2,19,['y = ' num2str(b(2),3) ' radio + ' num2str(b(1),3)])

% newspaper
X = [ones(length(sales),1) newspaper];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

subplot 313
hold on
plot(newspaper,salesPredict,'r','linewidth',2)
hold off

R2 = stats(1);
Fstat = stats(2);
pValue = stats(3);
var = stats(4);

text(80,28,['R^2 = ' num2str(R2)])
text(80,22,['y = ' num2str(b(2),3) ' newspaper + ' num2str(b(1),3)])
%% two predictors

% TV - radio
X = [ones(length(sales),1) TV radio];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

figure
% hold on
plot3(TV,radio,sales,'.r','linewidth',2)
xlabel('TV')
ylabel('radio')
zlabel('sales')
grid on

% grid config
nPointsx=20;
stepx = (max(TV)-min(TV))/nPointsx;
xgrid=min(TV):stepx:max(TV);
nPointsy=20;
stepy = (max(radio)-min(radio))/nPointsy;
xgrid=min(TV):stepx:max(TV);
ygrid=min(radio):stepy:max(radio);

[X_grid, Y_grid] = meshgrid(xgrid,ygrid);

for ii = 1:length(xgrid)
    for jj = 1: length(ygrid)
        Z_grid(ii,jj) = [1 xgrid(ii) ygrid(jj)]*b;
    end
end

hold on
m = mesh(X_grid,Y_grid,Z_grid,'FaceAlpha',0.5);
hold off

R2 = stats(1);
title(['y = ' num2str(b(3),3) ' radio + ' num2str(b(2),3) ' TV + ' num2str(b(1),3) ' , ' 'R^2 = ' num2str(R2)])

%%
% xx=[];
% yy=[];
% for ii = 1:length(xgrid)
%     for jj = 1:length(ygrid)
%         xx=[xx X_grid(ii,jj)+Y_grid(ii,jj)];
%         yy =[yy Z_grid(ii,jj)];
%     end
% end
% 
% plot(xx,yy)
% 
% pause

% TV - newspaper
X = [ones(length(sales),1) TV newspaper];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

figure
% hold on
plot3(TV,newspaper,sales,'.r','linewidth',2)
xlabel('TV')
ylabel('newspaper')
zlabel('sales')
grid on

% grid config
nPointsx=20;
stepx = (max(TV)-min(TV))/nPointsx;
xgrid=min(TV):stepx:max(TV);
nPointsy=20;
stepy = (max(newspaper)-min(newspaper))/nPointsy;
xgrid=min(TV):stepx:max(TV);
ygrid=min(newspaper):stepy:max(newspaper);

[X_grid, Y_grid] = meshgrid(xgrid,ygrid);

for ii = 1:length(xgrid)
    for jj = 1: length(ygrid)
        Z_grid(ii,jj) = [1 xgrid(ii) ygrid(jj)]*b;
    end
end

hold on
m = mesh(X_grid,Y_grid,Z_grid,'FaceAlpha',0.5);
hold off

R2 = stats(1);
title(['y = ' num2str(b(3),3) ' newspaper + ' num2str(b(2),3) ' TV + ' num2str(b(1),3) ' , ' 'R^2 = ' num2str(R2)])


% radio - newspaper
X = [ones(length(sales),1) radio newspaper];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

figure
% hold on
plot3(radio,newspaper,sales,'.r','linewidth',2)
xlabel('radio')
ylabel('newspaper')
zlabel('sales')
grid on

% grid config
nPointsx=20;
stepx = (max(radio)-min(radio))/nPointsx;
xgrid=min(radio):stepx:max(radio);
nPointsy=20;
stepy = (max(newspaper)-min(newspaper))/nPointsy;
xgrid=min(radio):stepx:max(radio);
ygrid=min(newspaper):stepy:max(newspaper);

[X_grid, Y_grid] = meshgrid(xgrid,ygrid);

for ii = 1:length(xgrid)
    for jj = 1: length(ygrid)
        Z_grid(ii,jj) = [1 xgrid(ii) ygrid(jj)]*b;
    end
end

hold on
m = mesh(X_grid,Y_grid,Z_grid,'FaceAlpha',0.5);
hold off

R2 = stats(1);
title(['y = ' num2str(b(3),3) ' newspaper + ' num2str(b(2),3) ' radio + ' num2str(b(1),3) ' , ' 'R^2 = ' num2str(R2)])

%% interaction

%% two predictors

% TV - radio
X = [ones(length(sales),1) TV radio TV.*radio];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

figure
% hold on
plot3(TV,radio,sales,'.r','linewidth',2)
xlabel('TV')
ylabel('radio')
zlabel('sales')
grid on

% grid config
nPointsx=20;
stepx = (max(TV)-min(TV))/nPointsx;
xgrid=min(TV):stepx:max(TV);
nPointsy=20;
stepy = (max(radio)-min(radio))/nPointsy;
xgrid=min(TV):stepx:max(TV);
ygrid=min(radio):stepy:max(radio);

[X_grid, Y_grid] = meshgrid(xgrid,ygrid);

for ii = 1:length(xgrid)
    for jj = 1: length(ygrid)
        Z_grid(ii,jj) = [1 xgrid(ii) ygrid(jj) xgrid(ii)*ygrid(jj)]*b;
    end
end

hold on
m = mesh(X_grid,Y_grid,Z_grid,'FaceAlpha',0.5);
hold off

R2 = stats(1);
title(['y = ' num2str(b(4),3) ' radio + ' num2str(b(3),3) ' TV + ' num2str(b(2),3) ' TV x radio ' num2str(b(1),3) ' , ' 'R^2 = ' num2str(R2)])

%% three predictors

% radio - newspaper
X = [ones(length(sales),1) TV radio newspaper];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

R2 = stats(1);

%% Including outliers

figure
nn = [151 180];
sales_out = sales;
sales_out(nn) = sales_out(nn)+[1 -1]'*50;
% TV
X = [ones(length(sales),1) TV];
[b,bint,r,rint,stats] = regress(sales,X);
salesPredict = X*b;

plot(TV,sales,'o','markersize',4)

hold on
plot(TV,sales_out,'+','markersize',6)
xlabel('TV + outliers')
ylabel('sales')

plot(TV,salesPredict,'r','linewidth',2)
R2 = stats(1);
Fstat = stats(2);
pValue = stats(3);
var = stats(4);

text(10,25,['R^2 = ' num2str(R2)])
text(10,19,['y = ' num2str(b(2),3) ' TV + ' num2str(b(1),3)])

[b,bint,r,rint,stats] = regress(sales_out,X);
salesPredict_out = X*b;
plot(TV,salesPredict_out,'b','linewidth',2)
hold off

R2 = stats(1);
Fstat = stats(2);
pValue = stats(3);
var = stats(4);

text(10,25-25,['R^2 (including outliers) = ' num2str(R2)])
text(10,19-25,['y = ' num2str(b(2),3) ' TV + ' num2str(b(1),3)])