Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.
|
|
Math Forum
»
Discussions
»
sci.math.*
»
sci.math
Notice: We are no longer accepting new posts, but the forums will continue to be readable.
Topic:
Computing tables for use with Peirce's Criterion - in 1855 and 2008
Replies:
3
Last Post:
Nov 21, 2012 9:28 AM
|
 |
|
|
Re: Computing tables for use with Peirce's Criterion - in 1855 and 2008
Posted:
Nov 21, 2012 7:48 AM
|
|
Op zondag 5 oktober 2008 15:41:37 UTC+2 schreef Knud Thomsen het volgende: > In 1852 Peirce described his criterion for identifying outliers in a > sample 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 a > description of the algorithm he had used to calculate them (2). The > calculations were done manually, aided only by log-tables and some > special function tables he managed to get his hands on! > > Gould's first two tables list the (squares of the) critical deviation > x as a function of N, the number of values in the sample, and k, the > number of outliers being tested for. These tables differ with respect > the 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 k > psi x erfc(x/sqrt(2)) > -------------------------------- > > Next, two equations are formulated from which the critical deviation > can 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^N > > Gould's equation A' can be combined with the unlabeled one just above > equation 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)<Eps > > Although this works quite well on a computer (with some extra > termination logistics), much faster convergence can be achieved > through root-finding by Newton's Method. > I currently use the algorithm below (suitable for double precision > floating point). Some loop-invariant code have been kept explicit for > the sake of readability. > > ---------------------------------------------------------------------- > function PeirceDeviation(N,k,m) > // Peirce's critical deviation for N samples, k potential outliers and > m > // 'unknown quantities' (usually 1). > > if N-m-k<=0 then return(0) // result undefined > > LnQN = k*Ln(k) + (N-k)*Ln(N-k) - N*Ln(N) > > x = 1 > > repeat > 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-16 > > return(x) > ---------------------------------------------------------------------- > > Gould did not list x-values below 1, and when glancing down the tables > along the edge of the filled-in cells, one is tempted to ask if > something mathematically interesting can be said about the special > case 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..161P > http://adsabs.harvard.edu/abs/1852AJ......2..176P (erratum) > > 2. Gould, B.A.: On Peirce's Criterion for the Rejection of Doubtful > Observations, > with tables for facilitating its application. > Astronomical Journal, vol. 4, iss. 83, p. 81-87 (1855). > http://adsabs.harvard.edu/abs/1855AJ......4...81G > > 3. Acton, F.S.: Real Computing Made Real: Preventing Errors in > Scientific and > Engineering Calculations. Princeton University Press (1996). > > > Best regards, > > Knud Thomsen > Geologist, Denmark
Dear mr. Thomsen,
Thanks for your work. Sadly, I was unable to implement your pseudocode I posted a question about this on the mathematics stackExchange. Could you help me out?
http://math.stackexchange.com/questions/241958/how-to-compute-peirces-critical-deviations
|
|
|
|