% LINREG_EX  Linear regression example. It contains simple and multiple regression
% implementation and calculates regression coefficients, its confidence intervals
% and correlation matrices. It also computes R2, TSS, RSS, ESS and F. The script
% is compatible with GNU Octave and MATLAB.
%
% This script randomly generates the input data set and saves it in .csv file
% (either linreg2.csv or linreg3.csv)
%
% The output results of this script have been checked in StatGraphics Plus 5.0

% (C) 2015 Alexey Voskov
% License: public domain
function linreg_ex
%regtype = 'simple'; npoints = 20;
regtype = 'multiple'; npoints = 100;
if strcmp(regtype, 'simple')
  x = 3*rand(npoints, 1);
  y = 2 + 5*x + randn(size(x));
  func = @(b,x) (b(1) + b(2)*x);
  csvwrite('linreg2.csv', [x y]);
  xv = [0;1;2;3];
elseif strcmp(regtype, 'multiple')
  x1 = 2*rand(npoints, 1);
  x2 = 2*rand(size(x1));  
  x = [x1 x2];
  y = 3 + 4*x1 + 5*x2 + randn(size(x1));
  func = @(b,x) (b(1) + b(2)*x(:,1) + b(3)*x(:,2));
  X = [ones(size(x1)) x1 x2];
  csvwrite('linreg3.csv', [x1 x2 y]);
  xv = [0 0; 0.5 0.6; 1.4 1.5];
else
  error('Unknown regtype option');
end

% Calculate coefficients confidence intervals
X = [ones(size(y)) x];
%beta = inv(X'*X)*X'*y; % Variant 1
beta = X\y; %Variant 2
err = func(beta,x) - y;
sigma2 = err'*err / (numel(y)-numel(beta));
C = sigma2*inv(X'*X); % Covariance matrix!
sbeta = sqrt(diag(C));
dbeta = sbeta * tinv(0.975, numel(x) - numel(beta));
tbeta = beta ./ sbeta;
% Calculate correlation matrix for coefficients
nmat = repmat(sbeta(:),1,numel(beta));
rmat = C ./ (nmat .* nmat');
% Calculate some another coefficients
TSS = sum((y - mean(y)).^2);
RSS = err'*err;
ESS = TSS - RSS;
f_ess = numel(beta) - 1;
f_rss = numel(y) - numel(beta);
Femp = (ESS./f_ess) ./ (RSS./f_rss);
R2 = ESS./TSS;

% Type the information about the regresion results
fprintf('Linear regression results\n');
fprintf('%10s %10s %10s %7s\n', 'beta', 'sbeta', 'dbeta', 't');
fprintf('%10.5f %10.5f %10.5f %7.3f\n', [beta sbeta dbeta tbeta]');
fprintf('Standard error of estimation sigma: %.5f\n', sqrt(sigma2));
fprintf('Coefficient of determination R2:    %.5f\n', R2);
fprintf('Resudual sum of squares  (RSS): %.2f\n', RSS);
fprintf('Explained sum of squares (ESS): %.2f\n', ESS);
fprintf('Total sum of squares     (TSS): %.2f\n', TSS);
fprintf('Femp = %.3f; Fcrit = %.3f\n', Femp, finv(0.95, f_ess, f_rss));
fprintf('Correlation matrix for coefficient estimates\n');
matfmt = ['    ', repmat('%8.4f', 1, numel(beta)), '\n'];
fprintf(matfmt, rmat); % rmat is symmetric (no transposition is required)

% Calculate confidence and prediction intervals
[dypred, dyconf] = calc_dy(xv, sigma2, C, numel(y)-numel(beta));
ypred = func(beta,xv);
fprintf('%20s %10s %21s %21s\n', '', '', 'Prediction limits', 'Confidence limits');
fprintf('%20s %10s %10s %10s %10s %10s\n', ...
  'x', 'ypred', 'y-dy', 'y+dy', 'y-dy', 'y+dy');
for i = 1:size(xv,1)
  fprintf('%20s %10.5f %10.5f %10.5f %10.5f %10.5f\n', ...
    sprintf('%8.4f', xv(i,:)), ypred(i),...
      ypred(i) - dypred(i), ypred(i) + dypred(i), ...
      ypred(i) - dyconf(i), ypred(i) + dyconf(i));
end

% Show the result of modeling
close all;
if strcmp(regtype, 'simple')
  % a) Simple linear regression
  xv = (0:0.01:3)';
  ypred = func(beta,xv);
  [dypred, dyconf] = calc_dy(xv, sigma2, C, numel(y)-numel(beta));
  close all;
  hold on;
  plot(x,y,'bo','LineWidth',2);
  plot(xv, ypred, 'k-', 'LineWidth',2);
  plot(xv, ypred - dypred, 'b-', 'LineWidth', 2);
  plot(xv, ypred - dyconf, 'r-', 'LineWidth', 2);
  plot(xv, ypred + dypred, 'b-', 'LineWidth', 2);
  plot(xv, ypred + dyconf, 'r-', 'LineWidth', 2);
  set(gca, 'FontSize', 20);
  xlabel('x', 'FontSize', 24);
  ylabel('y', 'FontSIze', 24);
  legend('Data', 'Model', 'Pred.Int.', 'Conf.Int','Location','NorthWest');
  hold off;
elseif strcmp(regtype, 'multiple');
  % b) Multiple linear regression
  [x1m, x2m] = meshgrid(0:0.1:2, 0:0.1:2);
  ym = beta(1) + beta(2)*x1m + beta(3)*x2m;
  [dypred, dyconf] = calc_dy([x1m(:) x2m(:)], sigma2, C, numel(y)-numel(beta));
  dypred = reshape(dypred, size(x1m));
  dyconf = reshape(dyconf, size(x2m));
  plot3(x1, x2, y, 'bo', 'LineWidth', 2);
  hold on;
  surf(x1m, x2m, ym);
  dy = dyconf; % Check: predicted or confidence intervals on graph
  mesh(x1m, x2m, ym - dy, 'LineWidth', 2, 'AlphaData', 0.1);
  mesh(x1m, x2m, ym + dy, 'LineWidth', 2);
  set(gca, 'FontSize', 20);
  xlabel('x_1', 'FontSize', 24);
  ylabel('x_2', 'FontSize', 24);  
  zlabel('y', 'FontSIze', 24);
  box on;
  grid on;
  hold off;
end
print('linreg','-dpng', '-r150');
end

function [dypred, dyconf] = calc_dy(xv, sigma2, C, f)
% Calculate confidence and prediction intervals for given y values
% Inputs:
%   xv -- input x values
%   sigma2 -- regression standard error
%   C -- convariance matrix for regression coefficients
% Outputs:
%   dypred -- prediction intervals
%   dyconf -- confidence intervals

[m n] = size(xv);
dypred = nan(size(xv,1),1);
dyconf = nan(size(xv,1),1);
for i = 1:numel(dypred)
  x0 = [1;xv(i,:)'];
  dyconf2 = x0'*C*x0;
  dypred(i) = sqrt(sigma2 + dyconf2);
  dyconf(i) = sqrt(dyconf2);
end
dypred = dypred * tinv(0.975, f);
dyconf = dyconf * tinv(0.975, f);
end