% =======================================================
% Linear Regression example
% =======================================================
% Script: Victor M. Espinoza, vespinoza@uchile.cl
% =======================================================
% Santiago, Chile, 2018
% =======================================================

% close all
clc
clear

e = randn(100,1)*100;
x = (0:99)';
m = 0.2;
c = 1;
y0 = m*x.^2+c*x;

y = y0+e;

plot(x,y,'.','markersize',12)
xlabel('X')
ylabel('y')
title('Data')
box off
grid on
legend('Data','location','northeastoutside')

%% Estimating "b" coefficients using regress function (not the only one)
X = [ones(size(x)) x x.^2];
[b,bint,r,rint,stats] = regress(y,X);
yPredict = b(3)*x.^2+b(2)*x+b(1);

b

%% Ploting estimated model
hold on
plot(x,yPredict,...
    'r',...
    'linewidth',2)
xlabel('X')
ylabel('Y')
legend('Data','Estimated Model','location','northeastoutside')
box off

%% Ploting True model
plot(x,y0,...
    'color',[1 1 1]*0.55,...
    'linewidth',2)
hold off

legend('Data','Estimated Model','True Model','location','northeastoutside')

%% Ploting equations and stats
R2 = stats(1);
Fstat = stats(2);
pValue = stats(3);
var = stats(4);

text(5,1500,['R^2 = ' num2str(R2)])
text(5,2100,['y = ' num2str(b(3),3) 'x^2 + ' num2str(b(2),3) 'x + ' num2str(b(1),3)])

%% Other stats tools

figure(2)

% ploting residuals
subplot 221
scatter(x,r)
title('Residue')
xlabel('X')
ylabel('error')

% ploting estimations of p.d.f.'s
subplot 222
[counts, centers] = hist(r);
bar(centers, counts/sum(counts),'c')
[f, xi]=ksdensity(r);
hold on
scale = 7;
plot(xi,scale*f/sum(f),...
    'r',...
    'linewidth',2)
[par]=mle(r,'distribution','normal');
valuePdf = -4*par(2)+par(1):4*par(2)+par(1);
yPdf = pdf('normal',valuePdf,par(1),par(2));
plot(valuePdf,150*yPdf,'b','linewidth',2)
box off
hold off
xlabel('error')
ylabel('P.d.f.')
title('Histogram')

% Ploting boxplots
subplot 223
boxplot(r)
title('Boxplot')
xlabel('Residue')
ylabel('error')