[Rd] Kinderman-Ramage (PR#2846)

[Peter Dalgaard BSA]

leydold at statistik.wu-wien.ac.at writes:

> 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

Josef, I'm in the process of implementing your fixes, but while they
certainly improve things, I'm still seeing anomalies. Could you try
this (it takes a while, even on a fast machine)

m <- sapply(1:200,function(i)table(trunc(100*pnorm(rnorm(1e6)))))
plot(rowMeans(m), type="l")

and tell me if you see a double peak in the middle? It's only about
1.3% off-target, but given the accuracy of the constants in the
routine, it does seem a bit more than you should expect. I see the
same effect both with the Mersenne-Twister and with the Wichmann-Hill
uniform generators.

 
> 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)
> 
> ______________________________________________
> R-devel at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-devel
> 

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907