New default algorithm for R_unif_rand, used in sample, to resolve PR

17494.  Includes adding sample.kind to RNGkind to allow reverting to
the old algorithm for reproducing previous results. Thanks to Duncan
Murdoch for contributing the patch.


git-svn-id: https://svn.r-project.org/R/trunk@76160 00db46b3-68df-0310-9c12-caf00c1e9a41

Author: luke

doc/NEWS.Rd

@@ -23,6 +23,16 @@
 	 CMD build} are in serialization format version 2, and resave by
        default produces files in serialization format version 2 (unless
        the original is already in format version 3).
+
+       \item The default method for generating from a discrete uniform
+       distribution (used in \code{sample()}, for instance) has been
+       changed. This addresses the fact, pointed out by Ottoboni and
+       Stark, that the previous method made \code{sample()} noticeably
+       non-uniform on large populations.  See \PR{17494} for a
+       discussion. The previous method can be requested using
+       \code{RNGkind()} if necessary for reproduction of old results.
+       Thanks to Duncan Murdoch for contributing the patch and Gabe
+       Becker for further assistance.
      }
   }
 

src/include/R_ext/Random.h

@@ -1,6 +1,6 @@
 /*
  *  R : A Computer Language for Statistical Data Analysis
- *  Copyright (C) 1998-2016    The R Core Team
+ *  Copyright (C) 1998-2019    The R Core Team
  *
  *  This header file is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU Lesser General Public License as published by
@@ -53,6 +53,11 @@ typedef enum {
     KINDERMAN_RAMAGE
 } N01type;
 
+/* Different ways to generate discrete uniform samples */
+typedef enum {
+    ROUNDING,
+    REJECTION
+} Sampletype;
 
 void GetRNGstate(void);
 void PutRNGstate(void);

src/include/Rmath.h0.in

@@ -361,6 +361,7 @@ double R_pow_di(double, int);
 
 double	norm_rand(void);
 double	unif_rand(void);
+double  R_unif_index(double);
 double	exp_rand(void);
 #ifdef MATHLIB_STANDALONE
 void	set_seed(unsigned int, unsigned int);

src/library/base/R/RNG.R

@@ -1,7 +1,7 @@
 #  File src/library/base/R/RNG.R
 #  Part of the R package, https://www.R-project.org
 #
-#  Copyright (C) 1995-2014 The R Core Team
+#  Copyright (C) 1995-2019 The R Core Team
 #
 #  This program is free software; you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
@@ -21,14 +21,15 @@
 ## The available kinds are in
 ## ../../../include/Random.h  and ../../../main/RNG.c [RNG_Table]
 ##
-RNGkind <- function(kind = NULL, normal.kind = NULL)
+RNGkind <- function(kind = NULL, normal.kind = NULL, sample.kind = NULL)
 {
     kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper",
                "Mersenne-Twister", "Knuth-TAOCP", "user-supplied",
                "Knuth-TAOCP-2002", "L'Ecuyer-CMRG", "default")
     n.kinds <- c("Buggy Kinderman-Ramage", "Ahrens-Dieter", "Box-Muller",
                  "user-supplied", "Inversion", "Kinderman-Ramage",
 		 "default")
+    s.kinds <- c("Rounding", "Rejection", "default")
     do.set <- length(kind) > 0L
     if(do.set) {
 	if(!is.character(kind) || length(kind) > 1L)
@@ -51,19 +52,34 @@ RNGkind <- function(kind = NULL, normal.kind = NULL)
                     domain = NA)
          if(normal.kind == length(n.kinds) - 1L) normal.kind <- -1L
     }
-    r <- 1L + .Internal(RNGkind(i.knd, normal.kind))
-    r <- c(kinds[r[1L]], n.kinds[r[2L]])
-    if(do.set || !is.null(normal.kind)) invisible(r) else r
+
+    if(!is.null(sample.kind)) {
+        if(!is.character(sample.kind) || length(sample.kind) != 1L)
+            stop("'sample.kind' must be a character string of length 1")
+        sample.kind <- pmatch(sample.kind, s.kinds) - 1L
+        if(is.na(sample.kind))
+            stop(gettextf("'%s' is not a valid choice", sample.kind),
+                 domain = NA)
+        if (sample.kind == 0L)
+            warning("non-uniform 'Rounding' sampler used",
+                    domain = NA)
+         if(sample.kind == length(s.kinds) - 1L) sample.kind <- -1L
+    }
+    r <- 1L + .Internal(RNGkind(i.knd, normal.kind, sample.kind))
+    r <- c(kinds[r[1L]], n.kinds[r[2L]], s.kinds[r[3L]])
+    if(do.set || !is.null(normal.kind) || !is.null(sample.kind)) 
+	invisible(r) else r
 }
 
-set.seed <- function(seed, kind = NULL, normal.kind = NULL)
+set.seed <- function(seed, kind = NULL, normal.kind = NULL, sample.kind = NULL)
 {
     kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper",
                "Mersenne-Twister", "Knuth-TAOCP", "user-supplied",
                "Knuth-TAOCP-2002", "L'Ecuyer-CMRG", "default")
     n.kinds <- c("Buggy Kinderman-Ramage", "Ahrens-Dieter", "Box-Muller",
                  "user-supplied", "Inversion", "Kinderman-Ramage",
 		 "default")
+    s.kinds <- c("Rounding", "Rejection", "default")
     if(length(kind) ) {
 	if(!is.character(kind) || length(kind) > 1L)
 	    stop("'kind' must be a character string of length 1 (RNG to be used).")
@@ -85,6 +101,18 @@ set.seed <- function(seed, kind = NULL, normal.kind = NULL)
                  domain = NA)
          if(normal.kind == length(n.kinds) - 1L) normal.kind <- -1L
     }
+    if(!is.null(sample.kind)) {
+        if(!is.character(sample.kind) || length(sample.kind) != 1L)
+            stop("'sample.kind' must be a character string of length 1")
+        sample.kind <- pmatch(sample.kind, s.kinds) - 1L
+        if(is.na(sample.kind))
+            stop(gettextf("'%s' is not a valid choice", sample.kind),
+                 domain = NA)
+        if (sample.kind == 0L)
+            warning("non-uniform 'Rounding' sampler used",
+                    domain = NA)
+         if(sample.kind == length(s.kinds) - 1L) sample.kind <- -1L
+    }
     .Internal(set.seed(seed, i.knd, normal.kind))
 }
 
@@ -96,9 +124,11 @@ RNGversion <- function(vstr)
     if (length(vnum) < 2L)
 	stop("malformed version string")
     if (vnum[1L] == 0 && vnum[2L] < 99)
-        RNGkind("Wichmann-Hill", "Buggy Kinderman-Ramage")
+        RNGkind("Wichmann-Hill", "Buggy Kinderman-Ramage", "Rounding")
     else if (vnum[1L] == 0 || vnum[1L] == 1 && vnum[2L] <= 6)
-	RNGkind("Marsaglia-Multicarry", "Buggy Kinderman-Ramage")
+	RNGkind("Marsaglia-Multicarry", "Buggy Kinderman-Ramage", "Rounding")
+    else if (vnum[1L] <= 2 || vnum[1L] == 3 && vnum[2L] <= 5)
+	RNGkind("Mersenne-Twister", "Inversion", "Rounding")
     else
-	RNGkind("Mersenne-Twister", "Inversion")
+	RNGkind("Mersenne-Twister", "Inversion", "Rejection")
 }

src/library/base/man/Random.Rd

@@ -1,6 +1,6 @@
 % File src/library/base/man/Random.Rd
 % Part of the R package, https://www.R-project.org
-% Copyright 1995-2018 R Core Team
+% Copyright 1995-2019 R Core Team
 % Distributed under GPL 2 or later
 
 \name{Random}
@@ -27,9 +27,9 @@
 \usage{
 \special{.Random.seed <- c(rng.kind, n1, n2, \dots)}
 
-RNGkind(kind = NULL, normal.kind = NULL)
+RNGkind(kind = NULL, normal.kind = NULL, sample.kind = NULL)
 RNGversion(vstr)
-set.seed(seed, kind = NULL, normal.kind = NULL)
+set.seed(seed, kind = NULL, normal.kind = NULL, sample.kind = NULL)
 }
 \arguments{
   \item{kind}{character or \code{NULL}.  If \code{kind} is a character
@@ -39,6 +39,10 @@ set.seed(seed, kind = NULL, normal.kind = NULL)
   \item{normal.kind}{character string or \code{NULL}.  If it is a character
     string, set the method of Normal generation.  Use \code{"default"}
     to return to the \R default.  \code{NULL} makes no change.}
+  \item{sample.kind}{character string or \code{NULL}.  If it is a character
+    string, set the method of discrete uniform generation (used in 
+    \code{\link{sample}}, for instance).  Use \code{"default"} to return to 
+    the \R default.  \code{NULL} makes no change.}
   \item{seed}{a single value, interpreted as an integer, or \code{NULL}
     (see \sQuote{Details}).}
   \item{vstr}{a character string containing a version number,
@@ -153,6 +157,12 @@ set.seed(seed, kind = NULL, normal.kind = NULL)
   whenever it is selected (even if it is the current normal generator)
   and when \code{kind} is changed.
 
+  \code{sample.kind} can be \code{"Rounding"} or \code{"Rejection"},
+  or partial matches to these.  The former was the default in versions
+  prior to 3.6.0:  it made \code{\link{sample}} noticeably non-uniform
+  on large populations, and should only be used for reproduction of old
+  results.  See \PR{17494} for a discussion. 
+
   \code{set.seed} uses a single integer argument to set as many seeds
   as are required.  It is intended as a simple way to get quite different
   seeds by specifying small integer arguments, and also as a way to get
@@ -163,7 +173,8 @@ set.seed(seed, kind = NULL, normal.kind = NULL)
   called with \code{seed = NULL} it re-initializes (see \sQuote{Note})
   as if no seed had yet been set.
 
-  The use of \code{kind = NULL} or \code{normal.kind = NULL} in
+  The use of \code{kind = NULL}, \code{normal.kind = NULL} or 
+  \code{sample.kind = NULL} in
   \code{RNGkind} or \code{set.seed} selects the currently-used
   generator (including that used in the previous session if the
   workspace has been restored): if no generator has been used it selects
@@ -198,16 +209,17 @@ set.seed(seed, kind = NULL, normal.kind = NULL)
   element \emph{codes} the kind of RNG and normal generator. The lowest
   two decimal digits are in \code{0:(k-1)}
   where \code{k} is the number of available RNGs.  The hundreds
-  represent the type of normal generator (starting at \code{0}).
+  represent the type of normal generator (starting at \code{0}), and
+  the ten thousands represent the type of discrete uniform sampler.
 
   In the underlying C, \code{.Random.seed[-1]} is \code{unsigned};
   therefore in \R \code{.Random.seed[-1]} can be negative, due to
   the representation of an unsigned integer by a signed integer.
 
-  \code{RNGkind} returns a two-element character vector of the RNG and
-  normal kinds selected \emph{before} the call, invisibly if either
-  argument is not \code{NULL}.  A type starts a session as the default,
-  and is selected either by a call to \code{RNGkind} or by setting
+  \code{RNGkind} returns a three-element character vector of the RNG,
+  normal and sample kinds selected \emph{before} the call, invisibly if 
+  either argument is not \code{NULL}.  A type starts a session as the 
+  default, and is selected either by a call to \code{RNGkind} or by setting
   \code{.Random.seed} in the workspace.
 
   \code{RNGversion} returns the same information as \code{RNGkind} about
@@ -287,7 +299,8 @@ set.seed(seed, kind = NULL, normal.kind = NULL)
   \bold{34}, 198 and \bold{35}, 89.
   \doi{10.2307/2347988}.
 }
-\author{of RNGkind: Martin Maechler. Current implementation, B. D. Ripley}
+\author{of RNGkind: Martin Maechler. Current implementation, B. D. Ripley
+  with modifications by Duncan Murdoch.}
 
 \seealso{
   \code{\link{sample}} for random sampling with and without replacement.

src/library/datasets/man/crimtab.Rd

@@ -1,6 +1,6 @@
 % File src/library/datasets/man/crimtab.Rd
 % Part of the R package, https://www.R-project.org
-% Copyright 1995-2018 R Core Team
+% Copyright 1995-2019 R Core Team
 % Distributed under GPL 2 or later
 
 \name{crimtab}
@@ -121,7 +121,7 @@ zobs <- (h.mean - mean(d.hei[,"height"]))/h.sd
 
 # 6) Replace infinite values by +/- 6 as in Student's paper:
 
-zobs[infZ <- is.infinite(zobs)] # 3 of them
+zobs[infZ <- is.infinite(zobs)] # none of them 
 zobs[infZ] <- 6 * sign(zobs[infZ])
 
 # 7) Plot the distribution:

src/library/utils/R/example.R

@@ -1,7 +1,7 @@
 #  File src/library/utils/R/example.R
 #  Part of the R package, https://www.R-project.org
 #
-#  Copyright (C) 1995-2018 The R Core Team
+#  Copyright (C) 1995-2019 The R Core Team
 #
 #  This program is free software; you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
@@ -62,12 +62,12 @@ function(topic, package = NULL, lib.loc = NULL,
                     add = TRUE)
 	} else {
 	    oldRNG <- RNGkind()
-	    on.exit(RNGkind(oldRNG[1L], oldRNG[2L]), add = TRUE)
+	    on.exit(RNGkind(oldRNG[1L], oldRNG[2L], oldRNG[3L]), add = TRUE)
 	}
 	## set RNG
 	if(is.logical(setRNG)) { # i.e. == TRUE: use the same as R CMD check
 	    ## see share/R/examples-header.R
-	    RNGkind("default", "default")
+	    RNGkind("default", "default", "default")
 	    set.seed(1)
 	} else eval(setRNG)
     }

src/library/utils/man/example.Rd

@@ -1,6 +1,6 @@
 % File src/library/utils/man/example.Rd
 % Part of the R package, https://www.R-project.org
-% Copyright 1995-2015 R Core Team
+% Copyright 1995-2019 R Core Team
 % Distributed under GPL 2 or later
 
 \name{example}
@@ -44,7 +44,7 @@ example(topic, package = NULL, lib.loc = NULL,
     \code{setRNG = TRUE} sets the same state as
     \command{R CMD \link{check}} does for
     running a package's examples.  This is currently equivalent to
-    \code{setRNG = \{RNGkind("default", "default"); set.seed(1)\}}.}
+    \code{setRNG = \{RNGkind("default", "default", "default"); set.seed(1)\}}.}
   \item{ask}{logical (or \code{"default"}) indicating if
     \code{\link{devAskNewPage}(ask = TRUE)} should be called
     before graphical output happens from the example code.  The value