```Date: Oct 5, 2008 9:41 AM
Author: sales@kt-algorithms.com
Subject: Computing tables for use with Peirce's Criterion - in 1855 and 2008

In 1852 Peirce described his criterion for identifying outliers in asample from an otherwise normal distribution (1).In 1855 Gould, having being applying the criterion professionally,published a set of tables of critical deviations together with adescription of the algorithm he had used to calculate them (2). Thecalculations were done manually, aided only by log-tables and somespecial function tables he managed to get his hands on!Gould's first two tables list the (squares of the) critical deviationx as a function of N, the number of values in the sample, and k, thenumber of outliers being tested for. These tables differ with respectthe value for m, the 'number of unknown quantities', set to 1 and 2,respectively.Below I have tried to be consistent with Gould's notation.Here are the exceptions:Gould            Here--------------------------------Lower case n     kpsi x            erfc(x/sqrt(2))--------------------------------Next, two equations are formulated from which the critical deviationcan be found (x>0).Gould's equation D can be stated:  R = R1(x) = exp((x^2-1)/2) erfc(x/sqrt(2))Defining (Gould's equation B):  Q^N = k^k (N-k)^(N-k) / N^NGould's equation A' can be combined with the unlabeled one just aboveequation C into:  R = R2(x) = Q^(N/k) [(N-m-k x^2)/(N-m-k)]^(-(N-k)/(2k))Now a numerical solution for x can be found after equating R1 and R2.Gould's algorithm was based on the following pseudoconstant strategy(3), the initial value for R being my own suggestion:  R = 0.25  <Find x from: R2(x)=R>  repeat    R = R1(x)    oldx = x    <Find x from: R2(x)=R>  until abs(x-oldx)<EpsAlthough this works quite well on a computer (with some extratermination logistics), much faster convergence can be achievedthrough root-finding by Newton's Method.I currently use the algorithm below (suitable for double precisionfloating point). Some loop-invariant code have been kept explicit forthe sake of readability.----------------------------------------------------------------------function PeirceDeviation(N,k,m)// Peirce's critical deviation for N samples, k potential outliers andm// 'unknown quantities' (usually 1).if N-m-k<=0 then return(0)  // result undefinedLnQN = k*Ln(k) + (N-k)*Ln(N-k) - N*Ln(N)x = 1repeat  x = min(x, sqrt((N-m)/k) - 1e-10)  // R1(x) and R2(x)  R1 = exp((x*x-1)/2)*erfc(x/sqrt(2))  R2 = exp( (LnQN - 0.5*(N-k)*ln((N-m-k*x*x)/(N-m-k)))/k )  // R1'(x) and R2'(x)  R1d = x*R1 - sqrt(2/pi/exp(1))  R2d = x*(N-k)/(N-m-k*x*x)*R2  // Newton  oldx = x  x = oldx - (R1-R2)/(R1d-R2d)until abs(x-oldx) < N * 2E-16return(x)----------------------------------------------------------------------Gould did not list x-values below 1, and when glancing down the tablesalong the edge of the filled-in cells, one is tempted to ask ifsomething mathematically interesting can be said about the specialcase x=1.For that purpose we shall accept non-integer values for N and/or k.Setting R1(1)=R2(1) we get:  erfc(1/sqrt(2)) = Q^(N/k) = k [(N-k)^(N-k) / N^N]^(1/k)or, defining the ratio C = k/N:  erfc(1/sqrt(2)) = C (1-C)^(1/C-1)So, x=1 implies a fixed k/N ratio, which does not depend on m.C must be estimated numerically; here it is to 30 decimal places:  C = 0.58974 36745 10713 07248 06181 37574 ...which just happens to be rather close to the simple fraction:  23/39  = 0.58974 35897 43589 ...1. Peirce, B.: Criterion for the rejection of doubtful observations.Astronomical Journal, vol. 2, iss. 45, p. 161-163 (1852).http://adsabs.harvard.edu/abs/1852AJ......2..161Phttp://adsabs.harvard.edu/abs/1852AJ......2..176P     (erratum)2. Gould, B.A.: On Peirce's Criterion for the Rejection of DoubtfulObservations,with tables for facilitating its application.Astronomical Journal, vol. 4, iss. 83, p. 81-87 (1855).http://adsabs.harvard.edu/abs/1855AJ......4...81G3. Acton, F.S.: Real Computing Made Real: Preventing Errors inScientific andEngineering Calculations. Princeton University Press (1996).Best regards,Knud ThomsenGeologist, Denmark
```