% Regression example data pairs - nonlinear:
% 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_nonlinear_fits.xls';
sheet_number = 1;
data_range_1 = 'A6:A26';
data_range_2 = 'D6:D26';
x = xlsread(xls_filename,sheet_number,data_range_1);
y = xlsread(xls_filename,sheet_number,data_range_2);

x_squared = x.^2;
x_cubed = x.^3;

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

% Linear regression analysis:
% (choose only the x data column)
N = length(x);
X_linear = [ones(N,1) x];
b_linear = X_linear\y; % Excel Coefficients

% Linear fit
x_fit = [0:0.01:2]';
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:
% (choose both the x data column and the x_squared column)
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,:))));

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

% Cubic Regression Analysis:
% (choose the x, x_squared, and x_cubed columns)
X_cubic = [ones(length(x),1) x x_squared x_cubed];
b_cubic = X_cubic\y; % Excel Coefficients

% cubic fit for the plot:
y_fit_cubic = b_cubic(1) + b_cubic(2)*x_fit + ...
    b_cubic(3)*x_fit.^2 + b_cubic(4)*x_fit.^3;

% Excel ANOVA: SS (Sum of Squares)
SS_regression_cubic = b_cubic'*X_cubic'*y - 1/N*(ones(N,1)'*y)^2;
SS_total = y'*y - 1/N*(ones(N,1)'*y)^2;
SS_residual_cubic = SS_total - SS_regression_cubic;

% Excel Regression Statistics
Multiple_R_cubic = sqrt(SS_regression_cubic/SS_total);
R_Square_cubic = Multiple_R_cubic^2;
Adjusted_R_Square_cubic = 1-(1-Multiple_R_cubic^2)*(N-1)/...
    (N-1-(length(X_cubic(1,2:end))));
Standard_Error_cubic = sqrt(sum([b_cubic(1)+b_cubic(2)*x+b_cubic(3)*x.^2+b_cubic(4)*x.^3-y].^2)/...
    (N-length(X_cubic(1,:))));

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

% Output results
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));

fprintf('\n\nCubic Regression Analysis\n');
fprintf('Regression Statistics:\n');
fprintf('Multiple R \t\t\t%g\n',Multiple_R_cubic);
fprintf('R Square \t\t\t%g\n',R_Square_cubic);
fprintf('Adjusted R Square \t%g\n',Adjusted_R_Square_cubic);
fprintf('Standard Error \t\t%g\n',Standard_Error_cubic);
fprintf('Observations \t\t%g\n',Observations);

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

figure(1)
hold on;
plot(x,y,'ro') % raw data
plot(x_fit,y_fit_linear,'b-','LineWidth',2) % linear
plot(x_fit,y_fit_quadratic,'g-','LineWidth',2) % quadratic
plot(x_fit,y_fit_cubic,'r-','LineWidth',2) % cubic
plot(x,y,'ro','MarkerFaceColor','r')
legend('data pairs','linear fit','quadratic fit','cubic fit','Location','Best')
ylim([0 4.5])
xlabel('\itx','FontWeight','Bold')
ylabel('\ity','FontWeight','Bold')
grid on;
hold off;