%Example 3.13 in Holman: Chauvenet's Criterion

x=[5.30;5.73;6.77;5.26;4.33;5.45;6.09;5.64;5.81;5.75]; %sample points from Ex. 3.7

%Find and display sample size
n=length(x)

xm=mean(x);                                            %sample mean and estimate of population mean

sigm=std(x);                                           %Estimate of population standard deviation (uses n-1)

%Now calculate deviations from the mean, in absolute value:

devs=abs(x-xm);

%Find and display the ratio of the deviations to the std:

devs_ratio=devs/sigm

%Just compare to maximum allowable ratio according to the sample size
%(Table 3.5)

%If we wish to automate the process (hundreds of points), we can have Table
%3.5 here:

n_sample=[3 4 5 6 7 10 15 25 50 100 300 500 1000]';             %I used the transpose to make a column vector without typing the ;
max_dev_ratio=[1.38 1.54 1.65 1.73 1.80 1.96 2.13 2.33 2.57 2.81 3.14 3.29 3.48]';

%Now the detection routine (uses interpolation/extrapolation in case n is not listed)

maxdev=interp1(n_sample,max_dev_ratio,n);

k=0;
for i=1:n,
    if devs_ratio(i)>maxdev,
        i
        disp('Point eliminated')
    else
        k=k+1;
        remaining_x(k)=x(i);
    end
end

%Display remaining points

remaining_x