function [out] = walshnp(x, r, alpha)
%Checks for a nugget of r outliers, in the left or the right tail,
% where r is specified. For significance level  alpha = 0.05,
% the minimum sample size is 221.
% Example of use:
% 
%  > dat = normrnd(0, 2, [300 1]);
%  > data = [-11.26; dat; 13.12];
%  > walshnp(data)
%    %Lower outliers are: -11.26
%    %Upper outliers are: 13.12
%    %ans = -11.2600   13.1200

if nargin == 2
    alpha = 0.05
end
if nargin == 1
    alpha = 0.05; r = 1;
end
x = x(:);
x = sort(x);

n = length(x);
if n < 1/2 * (1 + 1/alpha)^2
    str = ['Impossible to achieve significance ',num2str(alpha),...
            ' with sample size n=', num2str(n)];
    error(str)
end

b2 = 1/alpha;
c = floor(sqrt(2 * n));
k = r + c;
a = (1 + sqrt( b2 * (c-b2)/(c-1)))/(c - b2 -1 );
rL = x(r) - (1 + a) * x(r+1) + a * x(k);
rU = x(n - r + 1) - (1 + a) * x(n-r) + a * x(n - k +1);

out = []; outL=[]; outU=[];
if (rL < 0) 
    outL = x(1:r); 
    str=['Lower outliers are: ', num2str(outL')];
    disp(str) 
    out = [out outL];
end
if (rU > 0) 
    outU = x(n-r+1:n);
    str=['Upper outliers are: ', num2str(outU')];
    disp(str) 
    out =[out outU];
end
if (rL > 0 & rU < 0)
    str =['No outliers when r = ',num2str(r)];
    disp(str)
end

