The R Project SVN R

Rev

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