[Rd] Kinderman-Ramage (PR#2846)

[leydold at statistik.wu-wien.ac.at]

Hi,

Our department has detected a bug in the implementation of the
Kinderman-Ramage generator for normal random variates in version
1.7.0, which can be seen from the below R session.
(Consecutive calls for chisq.test(...) always gives p-values very
close to 0.)
We have already encountered this bug in version 1.6.2

The error is in file 
R-1.7.0/src/nmath/snorm.c

Here is a patch for this file to correct the error
(the added two lines are crucial):

--- R-1.7.0/src/nmath/snorm.c	Wed Feb 26 16:51:17 2003
+++ snorm.c	Fri Apr 25 09:22:28 2003
@@ -197,7 +197,9 @@
 	u1 = unif_rand();
 	if(u1 < 0.884070402298758) {
 	    u2 = unif_rand();
-	    return A*(1.13113163544180*u1+u2-1);
+	    /** <modified constant> **/
+	    return A*(1.131131635444180*u1+u2-1);
+	    /** </modified> **/
 	}
 	
 	if(u1 >= 0.973310954173898) { /* tail: */
@@ -241,6 +243,10 @@
 	    tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
 	    if(fmax2(u2,u3) <= 0.805577924423817)
 		return (u2<u3) ? tt : -tt;
+	    /** <added> **/
+	    if(0.053377549506886*fabs(u2-u3) <= g(tt))
+	      return (u2<u3) ? tt : -tt;
+	    /** </added> **/
 	}
     case BOX_MULLER:
 	if(BM_norm_keep != 0.0) { /* An exact test is intentional */



----------------------------------------------------------

$ uname -a
Linux 2.4.18-3 #1 Thu Apr 18 07:37:53 EDT 2002 i686 unknown

$ R
R : Copyright 2003, The R Development Core Team
Version 1.7.0  (2003-04-16)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

> RNGkind(normal.kind="Kinderman-Ramage")
> chisq.test(table(trunc(100*pnorm(rnorm(1e6)))))

        Chi-squared test for given probabilities

data:  table(trunc(100 * pnorm(rnorm(1e+06)))) 
X-squared = 204.1378, df = 99, p-value = 2.746e-09


Josef

--

-----------------------------------------------------------------------------
Josef Leydold     |     University of Economics and Business Administration
                  |     Department for Applied Statistics and Data Processing
-----------------------------------------------------------------------------
Augasse 2-6       |     Tel.   *43/1/31336-4695
A-1090 Vienna     |     FAX    *43/1/31336-738
European Union    |     email  Josef.Leydold at statistik.wu-wien.ac.at
-----------------------------------------------------------------------------
Alles Unglueck kam daher, dass die Denkenden nicht mehr handeln konnten,
und die Handelnden keine Zeit mehr fanden zu denken.       (Marlen Haushofer)