| 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}
|