function ANOVAdemo
% ANOVAdemo  Example of a one-way ANOVA.  Uses data from D.C. Montomery
%            _Design and Analysis of Experiments_, 5th ed., 2001,
%            Wiley, New York.  See Chapter 3, pp. 60-62
%
%  No inputs to the main program

[Seq,CWP,TS] = readCSVData('Montgomery_Cotton_Data_set.csv',3,1);
ANOVA_oneway(TS,CWP)


% ====================================================
function ANOVA_oneway(x,levels,verbose)
% ANOVA_oneway  One-way ANOVA
%
% Synopsis:  ANOVA_oneway(x,levels)
%            ANOVA_oneway(x,levels,verbose)
%
% Input: x = vector of response data for multiple treatments of a single factor
if nargin<3, verbose=true;  end

% -- Find the unique levels
uL = unique(levels,'first');   %  Vector of unique levels
nL = length(uL);               %  Number of uniqe levels

% -- Extract data for each factor and compute ANOVA sums on the fly
xbar = zeros(1,nL);      %  mean for each treatment
nxL = zeros(1,nL);       %  Number of samples in each treatment
n = length(x);           %  Total number of samples (all levels)
Xbarbar = mean(x(:));    %  The grand mean
SSA = 0;   SSW = 0;   SST = 0;   %  Initialize sums to zero
for j=1:length(uL)
  xL = x( levels == uL(j) );    %  x values for this level
  nxL(j) = length( xL );
  xbar(j) = mean( xL );
  SSA = SSA + nxL(j)*(xbar(j) - Xbarbar)^2;
  SSW = SSW + sum( ( xL - xbar(j) ).^2 );
  SST = SST + sum( ( xL - Xbarbar ).^2 );
end

% -- Check for equal length data
if any(nxL~=nxL(1))
  warning('ANOVA:UnbalancedFtest','Unequal number of data for each factor');
  if verbose
    fprintf('\nSizes of samples for each factor\n   i          N\n');
    for j=1:length(nxL)
      fprintf('%4d  %9d\n',j,nxL(j));
    end
  end
end

% -- Compute summary statistics and test whether F > Fcrit
c = nL;              %  Use c to match formulas in Levine, Ramsey & Smidt
MSA = SSA/(c - 1);
MSW = SSW/(n - c);
F = MSA/MSW;                       %  F statistic for the test
alpha = 0.05;                      %  Significance level
fcrit = fq( 1-alpha, c-1, n-c );   %  Critical value of F, from Statbox toolbox
pval = 1 - fp( F, c-1, n-c );      %  Probability of F by chance; from Statbox

fprintf('\nANOVA summary\n');
fprintf('                          Sum of Sq.      Mean Square\n');
fprintf('Source           DOF       Squares        (Variance)      F\n%s\n',repmat('-',1,60));
fprintf('Among groups:  %4d  %14.5f  %14.5f  %8.3f\n',c-1,SSA,MSA,F)
fprintf('Within groups: %4d  %14.5f  %14.5f\n',n-c,SSW,MSW)
fprintf('Total:         %4d  %14.5f\n',n,SST);
fprintf('\nfcrit, pval = %12.4f  %12.3e\n',fcrit,pval);
fprintf('\nChecksum:  SST - (SSA+SSW) = %12.3e\n',SST-(SSA+SSW));

fprintf('\n Data for each treatment:\n Xbarbar = %14.5f\n\n',Xbarbar);
fprintf('   j    nxL      xbar(j)\n')
for j=1:nL
  fprintf(' %3d   %3d   %11.4f\n',j,nxL(j),xbar(j));
end

% ====================================================
function [x,y,z] = readCSVData(fname,ncol,nhead)
% readCSVData  Import data from a CSV file

if nargin<1, fname='testNumText.dat';  end
if nargin<2, ncol = 2;         end
if nargin<3, nhead = 0;        end  %  number of header lines
if nargin<4, verbose = false;  end

% --- Read entire file into cell array C
fid = fopen(fname);
fmt = repmat('%f',1,ncol);
C = textscan(fid, fmt,'delimiter',',','headerlines',nhead);
fclose(fid);

% --- Move data from elements of C to vectors.  C{i} is column i
x = C{1};
y = C{2};
z = C{3};

if verbose
  fprintf('\nData read from %s:\n',fname);
  for i=1:length(x)
    fprintf('%12.3e  %12.3e  %12.3e\n',x(i),y(i),z(i))
  end
end