% Example to be used in class to illustrate 
% a second-order polynomial curve fit
% Barukyah Shaparenko and John M. Cimbala

clear all
close all

% Import data from Excel
% NOTE: Excel data file must reside in same folder as Matlab m-file
% Use xlsread or xlsxread

xls_filename = 'Regression_second_order.xls';
sheet_number = 1;
data_range_1 = 'A5:A15';
data_range_2 = 'B5:B15';
data_range_3 = 'C5:C15';
x = xlsread(xls_filename,sheet_number,data_range_1);
x_squared = xlsread(xls_filename,sheet_number,data_range_2);
y = xlsread(xls_filename,sheet_number,data_range_3);

% -----------------------------------------------

% Linear Regression Analysis:
N = length(x);
X_linear = [ones(N,1) x];
b_linear = X_linear\y; % Excel Coefficients

% Straight line for the plot:
x_fit = [0:0.01:1]';
y_fit_linear = b_linear(1) + b_linear(2)*x_fit;

% Excel ANOVA: SS (Sum of Squares)
SS_regression_linear = b_linear'*X_linear'*y - 1/N*(ones(N,1)'*y)^2;
SS_total = y'*y - 1/N*(ones(N,1)'*y)^2;
SS_residual_linear = SS_total - SS_regression_linear;

% Excel Regression Statistics
Multiple_R_linear = sqrt(SS_regression_linear/SS_total);
R_Square_linear = Multiple_R_linear^2;
Adjusted_R_Square_linear = 1-(1-Multiple_R_linear^2)*(N-1)/(N-1-(length(X_linear(1,2:end))));
Standard_Error_linear = sqrt(sum([b_linear(1)+b_linear(2)*x-y].^2)/(N-length(X_linear(1,:))));
Observations = length(x);

% -----------------------------------------------

% Quadratic Regression Analysis:
X_quadratic = [ones(N,1) x x_squared];
b_quadratic = X_quadratic\y; % Excel Coefficients

% Quadratic fit for the plot:
y_fit_quadratic = b_quadratic(1) + b_quadratic(2)*x_fit + ...
    b_quadratic(3)*x_fit.^2;

% Excel ANOVA: SS (Sum of Squares)
SS_regression_quadratic = b_quadratic'*X_quadratic'*y - 1/N*(ones(N,1)'*y)^2;
SS_total = y'*y - 1/N*(ones(N,1)'*y)^2;
SS_residual_quadratic = SS_total - SS_regression_quadratic;

% Excel Regression Statistics
Multiple_R_quadratic = sqrt(SS_regression_quadratic/SS_total);
R_Square_quadratic = Multiple_R_quadratic^2;
Adjusted_R_Square_quadratic = 1-(1-Multiple_R_quadratic^2)*(N-1)/...
    (N-1-(length(X_quadratic(1,2:end))));
Standard_Error_quadratic = sqrt(sum([b_quadratic(1)+b_quadratic(2)*x+b_quadratic(3)*x.^2-y].^2)/...
    (N-length(X_quadratic(1,:))));

% -----------------------------------------------

fprintf('\n----------------------\n');
fprintf('Linear Regression Analysis\n');
fprintf('Regression Statistics:\n');
fprintf('Multiple R \t\t\t%g\n',Multiple_R_linear);
fprintf('R Square \t\t\t%g\n',R_Square_linear);
fprintf('Adjusted R Square \t%g\n',Adjusted_R_Square_linear);
fprintf('Standard Error \t\t%g\n',Standard_Error_linear);
fprintf('Observations \t\t%g\n',Observations);

fprintf('\nLinear Regression Analysis\n');
fprintf('Intercept = \t%g\n',b_linear(1));
fprintf('X Variable 1 = \t%g\n',b_linear(2));

fprintf('\n\nQuadratic Regression Analysis\n');
fprintf('Regression Statistics:\n');
fprintf('Multiple R \t\t\t%g\n',Multiple_R_quadratic);
fprintf('R Square \t\t\t%g\n',R_Square_quadratic);
fprintf('Adjusted R Square \t%g\n',Adjusted_R_Square_quadratic);
fprintf('Standard Error \t\t%g\n',Standard_Error_quadratic);
fprintf('Observations \t\t%g\n',Observations);

fprintf('\nQuadratic Regression Analysis\n');
fprintf('Intercept = \t%g\n',b_quadratic(1));
fprintf('X Variable 1 = \t%g\n',b_quadratic(2));
fprintf('X Variable 2 = \t%g\n',b_quadratic(3));

figure(1)
hold on;
plot(x,y,'ko') % raw data
plot(x_fit,y_fit_linear,'k--','LineWidth',2) % linear
plot(x_fit,y_fit_quadratic,'r-','LineWidth',2) % quadratic
legend('raw data','linear','quadratic','Location','South','Orientation','Horizontal')
plot(x,y,'ko','MarkerFaceColor','w')
ylim([1 3])
xlabel('\itx','FontWeight','Bold')
ylabel('\ity','FontWeight','Bold')
grid on;
set(gca,'Xcolor','g','Ycolor','g')
c_axes = copyobj(gca,gcf);
set(c_axes,'color','none','xcolor','k','xgrid','off','ycolor','k','ygrid','off')
hold off;