From: "Charles C. Kankelborg" <kankel@solar.physics.montana.edu>
Date: Tue May 6, 2003  3:31:57  AM Canada/Mountain
To: bugs@nr.com
Cc: Charles Kankelborg <kankel@mithra.physics.montana.edu>
Subject: ks2d1s.c bug (v2.10)

I have found a bug affecting ks2d1s.c and ks2d2s.c (NR version 2.10). The routines' results do not match the Monte-Carlo simulations published by Fasano & Franceschini (1987 MNRAS 225, 155). The essence of the problem is that the NR code assigns the present point (x,y) arbitrarily to a particular quadrant. The proper behavior (verified in the example below) is to choose the quadrant that maximizes the D-statistic. This is just the 2D generalization of what is done when comparing cumulative distributions in the 1D Kolmogorov-Smirnov test: one should always calculate D immediately to the left and immediately to the right of x. The bug fix is accomplished with only a few added lines of code (see below).

Enclosed please find:

ks2d_test1.c --- test code to replicate published Monte-Carlo results.
quadct.c --- revised NR routine
ks2d1s.c --- revised NR routine

Similar corrections should be made to ks2d2s.c. When both revised routines are used, the Fasano & Franceschini results (specifically Table A1) are reproduced accurately. The difference is most significant for small n. Here are two comparisons, using uncorrelated uniform distributions in x and y:

n = 5
-----
Test of ks2d1s (sample vs. uniform distro), using NR version 2.10:
Of 100000 trials with 5 points, 30 % had Zn < 0.683075 (prob=0.698173)
Of 100000 trials with 5 points, 40 % had Zn < 0.743544 (prob=0.528061)
Of 100000 trials with 5 points, 50 % had Zn < 0.803929 (prob=0.435169)
Of 100000 trials with 5 points, 60 % had Zn < 0.867547 (prob=0.369958)
Of 100000 trials with 5 points, 70 % had Zn < 0.943754 (prob=0.238967)
Of 100000 trials with 5 points, 80 % had Zn < 1.037425 (prob=0.164640)
Of 100000 trials with 5 points, 90 % had Zn < 1.167895 (prob=0.080599)
Of 100000 trials with 5 points, 95 % had Zn < 1.284298 (prob=0.043063)
Of 100000 trials with 5 points, 99 % had Zn < 1.512351 (prob=0.008434)

Test of ks2d1s with CCK modified NR routines:
Of 100000 trials with 5 points, 30 % had Zn < 0.881128 (prob=0.313280)
Of 100000 trials with 5 points, 40 % had Zn < 0.941218 (prob=0.250265)
Of 100000 trials with 5 points, 50 % had Zn < 0.999598 (prob=0.191043)
Of 100000 trials with 5 points, 60 % had Zn < 1.059472 (prob=0.142157)
Of 100000 trials with 5 points, 70 % had Zn < 1.124804 (prob=0.098314)
Of 100000 trials with 5 points, 80 % had Zn < 1.205026 (prob=0.064659)
Of 100000 trials with 5 points, 90 % had Zn < 1.328328 (prob=0.029608)
Of 100000 trials with 5 points, 95 % had Zn < 1.440247 (prob=0.014036)
Of 100000 trials with 5 points, 99 % had Zn < 1.637874 (prob=0.006076)

n = 50
------
Test of ks2d1s (sample vs. uniform distro), using NR version 2.10:
Of 100000 trials with 50 points, 30 % had Zn < 0.904652 (prob=0.557122)
Of 100000 trials with 50 points, 40 % had Zn < 0.966805 (prob=0.472514)
Of 100000 trials with 50 points, 50 % had Zn < 1.028927 (prob=0.391126)
Of 100000 trials with 50 points, 60 % had Zn < 1.094240 (prob=0.317008)
Of 100000 trials with 50 points, 70 % had Zn < 1.169169 (prob=0.246135)
Of 100000 trials with 50 points, 80 % had Zn < 1.261041 (prob=0.175808)
Of 100000 trials with 50 points, 90 % had Zn < 1.398830 (prob=0.100489)
Of 100000 trials with 50 points, 95 % had Zn < 1.517613 (prob=0.058069)
Of 100000 trials with 50 points, 99 % had Zn < 1.757962 (prob=0.017663)

Test of ks2d1s with CCK modified NR routines (agrees with F&F 1987):
Of 100000 trials with 50 points, 30 % had Zn < 0.972676 (prob=0.464874)
Of 100000 trials with 50 points, 40 % had Zn < 1.034577 (prob=0.384917)
Of 100000 trials with 50 points, 50 % had Zn < 1.095263 (prob=0.317005)
Of 100000 trials with 50 points, 60 % had Zn < 1.161283 (prob=0.254064)
Of 100000 trials with 50 points, 70 % had Zn < 1.235185 (prob=0.193772)
Of 100000 trials with 50 points, 80 % had Zn < 1.326690 (prob=0.135520)
Of 100000 trials with 50 points, 90 % had Zn < 1.462602 (prob=0.075323)
Of 100000 trials with 50 points, 95 % had Zn < 1.581146 (prob=0.043305)
Of 100000 trials with 50 points, 99 % had Zn < 1.819722 (prob=0.011908)

In the above, Zn = D * sqrt(n) and prob is the value returned by the NR routine using the emprirical formula and Pearson's r calculated from the data. Note that since r is calculated from the data rather than from the user-supplied distribution function, the values of 1-prob are not expected to agree exactly with the confidence levels calculated by the Monte-Carlo.

Best regards, and keep up the good work!
Charles
________________________________________________________________________
Charles Kankelborg AC7NY |
260C EPS Building        | "Truth and falsehood are opposed; but truth
Physics Department       |  is the norm not of truth only but of
Montana State University |  falsehood also."
Bozeman, MT 59717        |
Phone: 406-994-7853      |          --C. S. Lewis, The Allegory of Love
FAX:   406-994-4452      |
http://solar.physics.montana.edu/kankel/
¨
¨
¨