| grubbs(x, x_date, sided, alpha, dist, nmin)
|
function [Gmax_loop Gmin_loop G2_loop lambda1_loop lambda2_loop outlier outlier_num] = grubbs(x, x_date, sided, alpha, dist, nmin)
%
% The 'grubbs' function performs an iterative test to check for the max,
% min or both observations as outliers. The test checks first for the max
% or min data point to be outlier and if any of the two observations is
% labeled as outlier then removes it from the sample and checks for the
% next max or min. If both of the test for max or min are not passed, then
% it checks if both max and min data points can be labeled as outliers. The
% test stops when no outliers are identified or the data sample is less
% than 'nmin'. The first part of the algorithm is also called the 'Modified
% Thompson Tau' procedure.
%
% Data in 'x' are organized so that columns are the time series and rows
% are the time intervals. All series contain the same number of
% observations.
%
% 'x_date' is a column vector (cell array) containing the dates of the
% corresponding data in 'x'.
%
% 'sided' indicates if the critical values are from a two-sided
% distribution (default) or one-sided ('sided' = 0).
%
% 'alpha' specifies the significant level. By default 'alpha' = 0.05.
%
% The variable 'dist' indicates if the series are assumed to be normal
% (default) or log-normally distributed. It can take the values 0 (normal)
% or 1 (log-normal).
%
% 'nmin' specifies the minimum number of observations, either in the
% original series or after removing outliers, before the algorithm stops.
% The default value is 6.
%
% [Gmax_loop Gmin_loop G2_loop lambda1_loop lambda2_loop outlier outlier_num outlier] =
% = gesd(...) returns the following:
% Gmax_loop - indicates the statistic for the max observation as
% outlier in each loop (rows).
% Gmin_loop - indicates the statistic for the min observation as
% outlier in each loop (rows).
% G2_loop - indicates the statistic for the max and min
% observations as outliers in each loop (rows).
% lambda1_loop - critical values for the max or mmin observation as
% outlier in each loop (rows).
% lambda2_loop - critical values for both the max and mmin
% observation as outliers in each loop (rows).
% outlier - cell array specifying the date and the series number (column
% number in 'x') where the potential outliers are situated.
% outlier_num - matrix providing row and column numbers of the values in
% 'x' considered as potential outliers.
%
% Created by Francisco Augusto Alcaraz Garcia
% alcaraz_garcia@yahoo.com
%
% References:
%
% 1) H.A David; H.O. Hartley; E.S. Pearson (1954). The Distribution of
% the Ratio, in a Single Normal Sample, of Range to Standard Deviation.
% Biometrika 41(3), pp. 482-493.
%
% 2) F.E. Grubbs (1969). Procedures for Detecting Outlying Observations
% in Samples. Technometrics 11(1), pp. 1-21.
%
% 3) NIST/SEMATECH e-Handbook of Statistical Methods (2010),
% http://www.itl.nist.gov/div898/handbook/.
%
% 4) S.L.R. Ellison; V.J. Barwick; T.J.D. Farrant (2009). Practical
% Statistics for the Analytical Scientist. A Bench Guide. RSC Publishing,
% 2nd ed., Cambridge.
% Check number of input arguements
if (nargin < 2) || (nargin > 6)
error('Requires two to six input arguments.')
end
% Define default values
if nargin == 2,
sided = 1;
alpha = 0.05;
dist = 0;
nmin = 6;
elseif nargin == 3,
alpha = 0.05;
dist = 0;
nmin = 6;
elseif nargin == 4,
dist = 0;
nmin = 6;
elseif nargin == 5,
nmin = 6;
end
% Normal transformation
if dist == 1,
x = log(x);
end
if sided == 1,
alpha1 = alpha/2;
end
% Check for validity of inputs
if ~isnumeric(x) || ~isreal(x) || ~iscellstr(x_date),
error('Input x must be a numeric array and x_date must be a string table.')
elseif alpha <= 0 || alpha >= 1,
error('The confidence level must be between zero and one.')
end
[n, c] = size(x);
xr = x;
j4 = 1;
j7 = 1;
j9 = 1;
outlier = [];
outlier_num = [];
Gmax_loop = [];
Gmin_loop = [];
G2_loop = [];
lambda1_loop = [];
lambda2_loop = [];
if n <= nmin,
error('Time series should contain more than %d data points.', nmin);
else
while (sum(j4) + sum(j7) + sum(j9))~= 0,
nr = n - sum(isnan(xr));
[i1,j1] = find(nr > nmin);
if isempty(j1),
disp(['Time series after removal of outliers contain less than ', num2str(nmin), ' data points.']);
break
end
p1 = 1 - alpha1./nr;
lambda1 = tinv(p1,nr-2).*(nr-1)./sqrt((nr-2+tinv(p1,nr-2).^2).*nr);
Gmax = (nanmax(xr) - nanmean(xr))./nanstd(xr);
[i2,j2] = find(xr == repmat(nanmax(xr),n,1)); % location of max
[i3,j3] = find(Gmax > lambda1);
j4 = ismember(j1,j3)';
Gmin = (nanmean(xr) - nanmin(xr))./nanstd(xr);
[i5,j5] = find(xr == repmat(nanmin(xr),n,1)); % location of min
[i6,j6] = find(Gmin > lambda1);
j7 = ismember(j1,j6)';
p2 = 1 - alpha./(nr.*(nr - 1));
lambda2 = sqrt((2.*(nr-1).*tinv(p2,nr-2).^2)./(n-2+tinv(p2,n-2).^2));
G2 = (nanmax(xr) - nanmin(xr))./nanstd(xr);
if sum(j4)~=0,
xr(i2(j4),j2(j4)) = NaN;
outlier = [outlier; x_date(i2(j4)) cellstr(strcat('Series', num2str(j2(j4))))];
outlier_num = [outlier_num; i2(j4) j2(j4)];
end
if sum(j7)~=0,
xr(i5(j7),j5(j7)) = NaN;
outlier = [outlier; x_date(i5(j7)) cellstr(strcat('Series', num2str(j5(j7))))];
outlier_num = [outlier_num; i5(j7) j5(j7)];
end
if sum(j4)==0 && sum(j7)==0,
[i8,j8] = find(G2 > lambda2);
j9 = ismember(j1,j8)';
if sum(j9)~=0,
xr(i2(j9),j2(j9)) = NaN;
xr(i5(j9),j5(j9)) = NaN;
outlier = [outlier; x_date(i2(j9)) cellstr(strcat('Series', num2str(j2(j9)))); ...
x_date(i5(j9)) cellstr(strcat('Series', num2str(j5(j9))))];
outlier_num = [outlier_num; i2(j9) j2(j9); i5(j9) j5(j9)];
end
j9 = 0;
end
Gmax_loop = [Gmax_loop; Gmax];
Gmin_loop = [Gmin_loop; Gmin];
G2_loop = [G2_loop; G2];
lambda1_loop = [lambda1_loop; lambda1];
lambda2_loop = [lambda2_loop; lambda2];
end
end
|
|
