[Rd] Kinderman-Ramage (PR#2846)
leydold at statistik.wu-wien.ac.at
leydold at statistik.wu-wien.ac.at
Fri Apr 25 14:19:04 MEST 2003
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)
More information about the R-devel
mailing list