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
%        levels = vector of labels for each level of treatment
%                 length(x) must equal length(levels)
%        verbose = (optional) flag to control amount of printed output

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 );      %  Number of samples for this level
  xbar(j) = mean( xL );       %  Mean for this treatment
  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