function [ind  outl grubbsstat crit] = grubbs(data, alpha)
% Usage [grubbsstat crit] = grubbs(data, alpha)
% Input: data - univariate data set
%            alpha- significance level
% Output:    grubbsstat - Grubbs statistics
%                   crit - critical value for level alpha  
% if grubbsstat > crit, the point generating the grubbsstat is an outlier
% if outlier is found, its index is ind and value outl=data(ind)
%if outliet is not found, index and value are both NaN
% Example of use:
% >> a=randn(1, 100); a(17)=6;
% >> [i au gr cr]=grubbs(a, 0.05)
%gives 
% i = 17
% au = 6
% gr = 5.3442
% cr = 3.3841
%--------------------------------------------------------------
%Grubbs test assumes Gaussian data


if nargin < 2 
    alpha = 0.05
end
grubbsstat = max(abs(data - mean(data)))/std(data);

n=length(data);
crit = (n-1)/sqrt(n) * sqrt( tinv(alpha/(2*n),n-2)^2/ ...
          (n-2 + tinv(alpha/(2*n), n-2)^2 ) );
     
      ind = NaN; outl=NaN;
if (grubbsstat > crit)
    absdiff = abs(data - mean(data));  
    ind = find(  max(abs(data - mean(data)))  ==  abs(data - mean(data)), 1 );
    outl = data(ind);
end