The R Project SVN R

Rev

Rev 76241 | Rev 80840 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
42333 ripley 1
% File src/library/base/man/Random.Rd
68948 ripley 2
% Part of the R package, https://www.R-project.org
76160 luke 3
% Copyright 1995-2019 R Core Team
42333 ripley 4
% Distributed under GPL 2 or later
5
 
2 r 6
\name{Random}
56186 murdoch 7
\alias{Random}
3076 pd 8
\alias{RNG}
9
\alias{RNGkind}
23377 pd 10
\alias{RNGversion}
7732 ripley 11
\alias{set.seed}
33493 ripley 12
\alias{.Random.seed}
15168 pd 13
\title{Random Number Generation}
2 r 14
\description{
7081 pd 15
  \code{.Random.seed} is an integer vector, containing the random number
25099 hornik 16
  generator (RNG) \bold{state} for random number generation in \R.  It
7799 maechler 17
  can be saved and restored, but should not be altered by the user.
3076 pd 18
 
19
  \code{RNGkind} is a more friendly interface to query or set the kind
7660 ripley 20
  of RNG in use.
21
 
23377 pd 22
  \code{RNGversion} can be used to set the random generators as they
25099 hornik 23
  were in an earlier \R version (for reproducibility).
46400 ripley 24
 
7790 ripley 25
  \code{set.seed} is the recommended way to specify seeds.
2 r 26
}
15168 pd 27
\usage{
54653 ripley 28
\special{.Random.seed <- c(rng.kind, n1, n2, \dots)}
15168 pd 29
 
76160 luke 30
RNGkind(kind = NULL, normal.kind = NULL, sample.kind = NULL)
23377 pd 31
RNGversion(vstr)
76160 luke 32
set.seed(seed, kind = NULL, normal.kind = NULL, sample.kind = NULL)
15168 pd 33
}
3076 pd 34
\arguments{
8204 ripley 35
  \item{kind}{character or \code{NULL}.  If \code{kind} is a character
46398 ripley 36
    string, set \R's RNG to the kind desired.  Use \code{"default"} to
37
    return to the \R default.  See \sQuote{Details} for the
47143 ripley 38
    interpretation of \code{NULL}.}
7799 maechler 39
  \item{normal.kind}{character string or \code{NULL}.  If it is a character
46398 ripley 40
    string, set the method of Normal generation.  Use \code{"default"}
41
    to return to the \R default.  \code{NULL} makes no change.}
76160 luke 42
  \item{sample.kind}{character string or \code{NULL}.  If it is a character
43
    string, set the method of discrete uniform generation (used in 
44
    \code{\link{sample}}, for instance).  Use \code{"default"} to return to 
45
    the \R default.  \code{NULL} makes no change.}
62311 ripley 46
  \item{seed}{a single value, interpreted as an integer, or \code{NULL}
47
    (see \sQuote{Details}).}
23377 pd 48
  \item{vstr}{a character string containing a version number,
76205 hornik 49
    e.g., \code{"1.6.2"}.  The default RNG configuration of the current
50
    \R version is used if \code{vstr} is greater than the current version.} 
7081 pd 51
  \item{rng.kind}{integer code in \code{0:k} for the above \code{kind}.}
46400 ripley 52
  \item{n1, n2, \dots}{integers.  See the details for how many are required
7081 pd 53
    (which depends on \code{rng.kind}).}
3076 pd 54
}
49649 ripley 55
%% source and more detailed references:  ../../../main/RNG.c
3076 pd 56
\details{
7799 maechler 57
  The currently available RNG kinds are given below.  \code{kind} is
58
  partially matched to this list.  The default is
23384 pd 59
  \code{"Mersenne-Twister"}.
7790 ripley 60
  \describe{
61
    \item{\code{"Wichmann-Hill"}}{
7081 pd 62
    The seed, \code{.Random.seed[-1] == r[1:3]} is an integer vector of
63
    length 3, where each \code{r[i]} is in \code{1:(p[i] - 1)}, where
64
    \code{p} is the length 3 vector of primes, \code{p = (30269, 30307,
65
      30323)}.
7790 ripley 66
    The Wichmann--Hill generator has a cycle length of
67
    \eqn{6.9536 \times 10^{12}}{6.9536e12} (=
68
    \code{prod(p-1)/4}, see \emph{Applied Statistics} (1984)
69
    \bold{33}, 123 which corrects the original article).}
7799 maechler 70
 
7790 ripley 71
    \item{\code{"Marsaglia-Multicarry"}:}{
7081 pd 72
    A \emph{multiply-with-carry} RNG is used, as recommended by George
7790 ripley 73
    Marsaglia in his post to the mailing list \file{sci.stat.math}.
7799 maechler 74
    It has a period of more than \eqn{2^{60}}{2^60} and has passed
8204 ripley 75
    all tests (according to Marsaglia).  The seed is two integers (all
7790 ripley 76
    values allowed).}
3076 pd 77
 
7790 ripley 78
    \item{\code{"Super-Duper"}:}{
3076 pd 79
    Marsaglia's famous Super-Duper from the 70's.  This is the original
80
    version which does \emph{not} pass the MTUPLE test of the Diehard
7081 pd 81
    battery.  It has a period of \eqn{\approx 4.6\times 10^{18}}{about
7799 maechler 82
      4.6*10^18} for most initial seeds.  The seed is two integers (all
7081 pd 83
    values allowed for the first seed: the second must be odd).
7799 maechler 84
 
65189 hornik 85
    We use the implementation by Reeds \emph{et al} (1982--84).
7799 maechler 86
 
7842 ripley 87
    The two seeds are the Tausworthe and congruence long integers,
7081 pd 88
    respectively.  A one-to-one mapping to S's \code{.Random.seed[1:12]}
89
    is possible but we will not publish one, not least as this generator
7790 ripley 90
    is \bold{not} exactly the same as that in recent versions of S-PLUS.}
7799 maechler 91
 
62069 hornik 92
    \item{\code{"Mersenne-Twister"}:}{
75569 luke 93
      From Matsumoto and Nishimura (1998); code updated in 2002.
94
      A twisted GFSR with period
18519 ripley 95
      \eqn{2^{19937} - 1}{2^19937 - 1} and equidistribution in 623
42961 ripley 96
      consecutive dimensions (over the whole period).  The \sQuote{seed} is a
18519 ripley 97
      624-dimensional set of 32-bit integers plus a current position in
98
      that set.
75569 luke 99
 
100
      R uses its own initialization method due to B. D. Ripley and is
101
      not affected by the initialization issue in the 1998 code of
102
      Matsumoto and Nishimura addressed in a 2002 update.
18519 ripley 103
    }
7799 maechler 104
 
62069 hornik 105
    \item{\code{"Knuth-TAOCP-2002"}:}{
46398 ripley 106
      A 32-bit integer GFSR using lagged Fibonacci sequences with
42261 ripley 107
      subtraction.  That is, the recurrence used is
18519 ripley 108
      \deqn{X_j = (X_{j-100} - X_{j-37}) \bmod 2^{30}%
109
      }{X[j] = (X[j-100] - X[j-37]) mod 2^30}
42961 ripley 110
      and the \sQuote{seed} is the set of the 100 last numbers (actually
18519 ripley 111
      recorded as 101 numbers, the last being a cyclic shift of the
112
      buffer).  The period is around \eqn{2^{129}}{2^129}.
42261 ripley 113
    }
46400 ripley 114
 
62069 hornik 115
    \item{\code{"Knuth-TAOCP"}:}{
42261 ripley 116
      An earlier version from Knuth (1997).
46400 ripley 117
 
42261 ripley 118
      The 2002 version was not backwards compatible with the earlier
119
      version: the initialization of the GFSR from the seed was altered.
120
      \R did not allow you to choose consecutive seeds, the reported
121
      \sQuote{weakness}, and already scrambled the seeds.
42547 ripley 122
 
123
      Initialization of this generator is done in interpreted \R code
124
      and so takes a short but noticeable time.
18519 ripley 125
    }
62069 hornik 126
    \item{\code{"L'Ecuyer-CMRG"}:}{
56570 ripley 127
      A \sQuote{combined multiple-recursive generator} from L'Ecuyer
128
      (1999), each element of which is a feedback multiplicative
129
      generator with three integer elements: thus the seed is a (signed)
130
      integer vector of length 6. The period is around
131
      \eqn{2^{191}}{2^191}.
18519 ripley 132
 
56570 ripley 133
      The 6 elements of the seed are internally regarded as 32-bit
134
      unsigned integers.  Neither the first three nor the last three
135
      should be all zero, and they are limited to less than
136
      \code{4294967087} and \code{4294944443} respectively.
61433 ripley 137
 
56570 ripley 138
      This is not particularly interesting of itself, but provides the
139
      basis for the multiple streams used in package \pkg{parallel}.
140
      % See \code{\link{RngStream}}.
141
    }
142
 
62069 hornik 143
    \item{\code{"user-supplied"}:}{
55555 ripley 144
      Use a user-supplied generator.  See \code{\link{Random.user}} for
42261 ripley 145
      details.
7842 ripley 146
    }
3076 pd 147
  }
7799 maechler 148
 
47402 ripley 149
  \code{normal.kind} can be \code{"Kinderman-Ramage"},
150
  \code{"Buggy Kinderman-Ramage"} (not for \code{set.seed}),
151
  \code{"Ahrens-Dieter"}, \code{"Box-Muller"}, \code{"Inversion"} (the
152
  default), or \code{"user-supplied"}.  (For inversion, see the
153
  reference in \code{\link{qnorm}}.)  The Kinderman-Ramage generator
66886 maechler 154
  used in versions prior to 1.7.0 (now called \code{"Buggy"}) had several
57176 ripley 155
  approximation errors and should only be used for reproduction of old
47402 ripley 156
  results.  The \code{"Box-Muller"} generator is stateful as pairs of
157
  normals are generated and returned sequentially.  The state is reset
158
  whenever it is selected (even if it is the current normal generator)
159
  and when \code{kind} is changed.
7799 maechler 160
 
76160 luke 161
  \code{sample.kind} can be \code{"Rounding"} or \code{"Rejection"},
162
  or partial matches to these.  The former was the default in versions
163
  prior to 3.6.0:  it made \code{\link{sample}} noticeably non-uniform
164
  on large populations, and should only be used for reproduction of old
165
  results.  See \PR{17494} for a discussion. 
166
 
62311 ripley 167
  \code{set.seed} uses a single integer argument to set as many seeds
7799 maechler 168
  as are required.  It is intended as a simple way to get quite different
7790 ripley 169
  seeds by specifying small integer arguments, and also as a way to get
170
  valid seed sets for the more complicated methods (especially
57176 ripley 171
  \code{"Mersenne-Twister"} and \code{"Knuth-TAOCP"}).  There is no
172
  guarantee that different values of \code{seed} will seed the RNG
62311 ripley 173
  differently, although any exceptions would be extremely rare.  If
174
  called with \code{seed = NULL} it re-initializes (see \sQuote{Note})
175
  as if no seed had yet been set.
46398 ripley 176
 
76160 luke 177
  The use of \code{kind = NULL}, \code{normal.kind = NULL} or 
178
  \code{sample.kind = NULL} in
46400 ripley 179
  \code{RNGkind} or \code{set.seed} selects the currently-used
46889 ripley 180
  generator (including that used in the previous session if the
181
  workspace has been restored): if no generator has been used it selects
182
  \code{"default"}.
3076 pd 183
}
8334 ripley 184
\note{
56607 ripley 185
  Initially, there is no seed; a new one is created from the current
65147 ripley 186
  time and the process ID when one is required.  Hence different
187
  sessions will give different simulation results, by default.  However,
188
  the seed might be restored from a previous session if a previously
189
  saved workspace is restored.
16060 maechler 190
 
8334 ripley 191
  \code{.Random.seed} saves the seed set for the uniform random-number
15168 pd 192
  generator, at least for the system generators.  It does not
193
  necessarily save the state of other generators, and in particular does
194
  not save the state of the Box--Muller normal generator.  If you want
46400 ripley 195
  to reproduce work later, call \code{set.seed} (preferably with
196
  explicit values for \code{kind} and \code{normal.kind}) rather than
197
  set \code{.Random.seed}.
25685 ripley 198
 
38402 ripley 199
  The object \code{.Random.seed} is only looked for in the user's
25685 ripley 200
  workspace.
46400 ripley 201
 
46351 ripley 202
  Do not rely on randomness of low-order bits from RNGs.  Most of the
203
  supplied uniform generators return 32-bit integer values that are
204
  converted to doubles, so they take at most \eqn{2^{32}}{2^32} distinct
46398 ripley 205
  values and long runs will return duplicated values (Wichmann-Hill is
46351 ripley 206
  the exception, and all give at least 30 varying bits.)
8334 ripley 207
}
2 r 208
\value{
7081 pd 209
  \code{.Random.seed} is an \code{\link{integer}} vector whose first
7842 ripley 210
  element \emph{codes} the kind of RNG and normal generator. The lowest
22051 ripley 211
  two decimal digits are in \code{0:(k-1)}
46398 ripley 212
  where \code{k} is the number of available RNGs.  The hundreds
76160 luke 213
  represent the type of normal generator (starting at \code{0}), and
214
  the ten thousands represent the type of discrete uniform sampler.
1032 maechler 215
 
7842 ripley 216
  In the underlying C, \code{.Random.seed[-1]} is \code{unsigned};
27625 ripley 217
  therefore in \R \code{.Random.seed[-1]} can be negative, due to
218
  the representation of an unsigned integer by a signed integer.
7799 maechler 219
 
76160 luke 220
  \code{RNGkind} returns a three-element character vector of the RNG,
221
  normal and sample kinds selected \emph{before} the call, invisibly if 
222
  either argument is not \code{NULL}.  A type starts a session as the 
223
  default, and is selected either by a call to \code{RNGkind} or by setting
76405 ripley 224
  \code{.Random.seed} in the workspace. (NB: prior to \R 3.6.0 the first
225
  two kinds were returned in a two-element character vector.)
7660 ripley 226
 
46398 ripley 227
  \code{RNGversion} returns the same information as \code{RNGkind} about
228
  the defaults in a specific \R version.
229
 
7660 ripley 230
  \code{set.seed} returns \code{NULL}, invisibly.
2 r 231
}
232
\references{
74265 hornik 233
  Ahrens, J. H. and Dieter, U. (1973).
42261 ripley 234
  Extensions of Forsythe's method for random sampling from the normal
235
  distribution.
74265 hornik 236
  \emph{Mathematics of Computation}, \bold{27}, 927--937.
42261 ripley 237
 
74265 hornik 238
  Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988).
24300 ripley 239
  \emph{The New S Language}.
74265 hornik 240
  Wadsworth & Brooks/Cole.
241
  (\code{set.seed}, storing in \code{.Random.seed}.)
24300 ripley 242
 
74265 hornik 243
  Box, G. E. P. and Muller, M. E. (1958).
244
  A note on the generation of normal random deviates.
245
  \emph{Annals of Mathematical Statistics}, \bold{29}, 610--611.
246
  \doi{10.1214/aoms/1177706645}.
2 r 247
 
74265 hornik 248
  De Matteis, A. and Pagnutti, S. (1993).
249
  Long-range Correlation Analysis of the Wichmann-Hill Random Number
250
  Generator.
251
  \emph{Statistics and Computing}, \bold{3}, 67--70.
252
  \doi{10.1007/BF00153065}.
7799 maechler 253
 
74265 hornik 254
  Kinderman, A. J. and Ramage, J. G. (1976).
42261 ripley 255
  Computer generation of normal random variables.
74265 hornik 256
  \emph{Journal of the American Statistical Association}, \bold{71},
257
  893--896.
258
  \doi{10.2307/2286857}.
42261 ripley 259
 
74265 hornik 260
  Knuth, D. E. (1997).
261
  \emph{The Art of Computer Programming}.
262
  Volume 2, third edition.\cr
42261 ripley 263
  Source code at \url{http://www-cs-faculty.stanford.edu/~knuth/taocp.html}.
264
 
74265 hornik 265
  Knuth, D. E. (2002).
266
  \emph{The Art of Computer Programming}.
267
  Volume 2, third edition, ninth printing.
61433 ripley 268
 
74265 hornik 269
  L'Ecuyer, P. (1999).
270
  Good parameters and implementations for combined multiple recursive
271
  random number generators.
272
  \emph{Operations Research}, \bold{47}, 159--164.
273
  \doi{10.1287/opre.47.1.159}.
42261 ripley 274
 
74265 hornik 275
  Marsaglia, G. (1997).
276
  \emph{A random number generator for C}.
277
  Discussion paper, posting on Usenet newsgroup \code{sci.stat.math} on
7790 ripley 278
  September 29, 1997.
7799 maechler 279
 
74265 hornik 280
  Marsaglia, G. and Zaman, A. (1994).
281
  Some portable very-long-period random number generators.
282
  \emph{Computers in Physics}, \bold{8}, 117--121.
283
  \doi{10.1063/1.168514}.
7790 ripley 284
 
74265 hornik 285
  Matsumoto, M. and Nishimura, T. (1998).
7799 maechler 286
  Mersenne Twister: A 623-dimensionally equidistributed uniform
7790 ripley 287
  pseudo-random number generator,
288
  \emph{ACM Transactions on Modeling and Computer Simulation},
289
  \bold{8}, 3--30.\cr
67395 ripley 290
  Source code formerly at \code{http://www.math.keio.ac.jp/~matumoto/emt.html}.\cr
291
  Now see \url{http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/c-lang.html}.
75013 maechler 292
 
74265 hornik 293
  Reeds, J., Hubert, S. and Abrahams, M. (1982--4).
294
  C implementation of SuperDuper, University of California at Berkeley.
295
  (Personal communication from Jim Reeds to Ross Ihaka.)
18519 ripley 296
 
74265 hornik 297
  Wichmann, B. A.  and Hill, I. D. (1982).
298
  Algorithm AS 183: An Efficient and Portable Pseudo-random Number
299
  Generator.
300
  \emph{Applied Statistics}, \bold{31}, 188--190; Remarks:
42261 ripley 301
  \bold{34}, 198 and \bold{35}, 89.
74265 hornik 302
  \doi{10.2307/2347988}.
2 r 303
}
76160 luke 304
\author{of RNGkind: Martin Maechler. Current implementation, B. D. Ripley
305
  with modifications by Duncan Murdoch.}
3076 pd 306
 
52774 ripley 307
\seealso{
308
  \code{\link{sample}} for random sampling with and without replacement.
61433 ripley 309
 
52774 ripley 310
  \link{Distributions} for functions for random-variate generation from
311
  standard distributions.
312
}
75013 maechler 313
\examples{
41508 ripley 314
require(stats)
315
 
75013 maechler 316
## Seed the current RNG, i.e., set the RNG status
317
set.seed(42); u1 <- runif(30)
318
set.seed(42); u2 <- runif(30) # the same because of identical RNG status:
319
stopifnot(identical(u1, u2))
2 r 320
 
75013 maechler 321
\donttest{## the default random seed is 626 integers, so only print a few
322
 runif(1); .Random.seed[1:6]; runif(1); .Random.seed[1:6]
323
 ## If there is no seed, a "random" new one is created:
324
 rm(.Random.seed); runif(1); .Random.seed[1:6]
325
}
43870 murdoch 326
ok <- RNGkind()
61168 ripley 327
RNGkind("Wich")  # (partial string matching on 'kind')
16060 maechler 328
 
25118 hornik 329
## This shows how 'runif(.)' works for Wichmann-Hill,
16060 maechler 330
## using only R functions:
331
 
2 r 332
p.WH <- c(30269, 30307, 30323)
8204 ripley 333
a.WH <- c(  171,   172,   170)
15168 pd 334
next.WHseed <- function(i.seed = .Random.seed[-1])
335
  { (a.WH * i.seed) \%\% p.WH }
2 r 336
my.runif1 <- function(i.seed = .Random.seed)
3076 pd 337
  { ns <- next.WHseed(i.seed[-1]); sum(ns / p.WH) \%\% 1 }
75013 maechler 338
set.seed(1998-12-04)# (when the next lines were added to the souRce)
2 r 339
rs <- .Random.seed
3076 pd 340
(WHs <- next.WHseed(rs[-1]))
341
u <- runif(1)
8893 maechler 342
stopifnot(
343
 next.WHseed(rs[-1]) == .Random.seed[-1],
344
 all.equal(u, my.runif1(rs))
345
)
3076 pd 346
 
347
## ----
348
.Random.seed
61168 ripley 349
RNGkind("Super") # matches  "Super-Duper"
3076 pd 350
RNGkind()
351
.Random.seed # new, corresponding to  Super-Duper
352
 
353
## Reset:
7766 ripley 354
RNGkind(ok[1])
30059 ripley 355
 
76241 maechler 356
RNGversion(getRversion()) # the default version for this R version
357
 
30059 ripley 358
## ----
46351 ripley 359
sum(duplicated(runif(1e6))) # around 110 for default generator
30059 ripley 360
## and we would expect about almost sure duplicates beyond about
61150 ripley 361
qbirthday(1 - 1e-6, classes = 2e9) # 235,000
75016 ripley 362
}
286 maechler 363
\keyword{distribution}
364
\keyword{sysdata}