% About simple and polynomial regression
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loading data and storing in vector y and x
%load income.txt
%y=income(:,2);
%x=income(:,1);
x=[150 55 220 130 95]';
y=[220 75 250 145 200]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%estimate the underlying coefficients
n=input('Enter the order of the polynormial: ');
X=[x.^0];
for i=1:n
    X=[X x.^i];
end
betahat=inv(X'*X)*X'*y%least squared solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculate some necessary statistics
yhat=X*betahat;
ymean=mean(y);
SSTO=(y-ymean*ones(length(y),1))'*(y-ymean*ones(length(y),1))
SSR=(yhat-ymean*ones(length(y),1))'*(yhat-ymean*ones(length(y),1))
SSE=(y-yhat)'*(y-yhat)
Rsquare=SSR/SSTO
subplot(2,2,1),plot(x,y,'m *')
xlabel('x')
ylabel('y')
title('Scatter plot : x v.s. y')
%axis([0,18,100,260])
hold on%to hold the current plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%make a plot of the regression line
xx=(min(x):0.5:max(x))';
XX=[xx.^0];
for i=1:n
    XX=[XX xx.^i];
end
    yy=XX*betahat;
plot(xx,yy,'r')
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%make residual plots:x vs. residuals
subplot(2,2,2),plot(x,y-yhat,'b *')
xlabel('x')
ylabel('residual')
title('Scatter plot : x v.s. residual')
%axis([0,18,-30,30])
hold on
x2=zeros(max(x),1);
plot(x2,'r -')
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%make residual plots:y_hat vs. residuals
subplot(2,2,3),plot(yhat,y-yhat,'g *')
xlabel('yhat')
ylabel('residual')
title('Scatter plot : yhat v.s. residual')
%axis([100,300,-30,30])
hold on
x3=zeros(ceil(max(yhat)),1);
plot(x3,'r -')
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% to make a prediction
check='N';
check=input('Do you want to make a prediction? Y:YES  N:NO (default)','s');
if check=='Y'
    x_value=input('Enter x value: ');
    y_value=polyval(flipud(betahat),x_value);
    fprintf('The predicted value for x=%8.2f is %8.2f',x_value,y_value);
end
