SIST ISO 28640:2010
Random variate generation methods
Random variate generation methods
This International Standard specifies methods for generating uniform and non-uniform random variates for Monte Carlo simulation purposes. Cryptographic random number generation methods are not included. This International Standard is applicable, inter alia, by - researchers, industrial engineers or experts in operations management, who use statistical simulation, - statisticians who need randomization related to SQC methods, statistical design of experiments or sample surveys, - applied mathematicians who plan complex optimization procedures that require the use of Monte Carlo methods, and - software engineers who implement algorithms for random variate generation.
Méthodes de génération de nombres pseudo-aléatoires
Metode generiranja naključnih spremenljivk
Ta mednarodni standard opredeljuje metode za generiranje enotnih in neenotnih naključnih spremenljivk za namene Monte Carlo simulacije. Kriptografske metode generiranja naključnih števil niso vključene. Ta mednarodni standard med drugim uporabljajo: – raziskovalci, industrijski inženirji ali strokovnjaki za upravljanje operacij, ki uporabljajo statistično simulacijo, - statistiki, ki potrebujejo doseganje naključnosti, povezano z metodami SQC, statističnim načrtovanjem eksperimentov ali vzorčnih raziskav, - uporabni matematiki, ki načrtujejo kompleksne postopke optimizacije, ki zahtevajo uporabo Monte Carlo metod,in - inženirji programske opreme, ki izvajajo algoritme za generiranje naključnih spremenljivk.
General Information
Standards Content (Sample)
INTERNATIONAL ISO
STANDARD 28640
First edition
2010-03-15
Random variate generation methods
Méthodes de génération de nombres pseudo-aléatoires
Reference number
ISO 28640:2010(E)
©
ISO 2010
---------------------- Page: 1 ----------------------
ISO 28640:2010(E)
PDF disclaimer
This PDF file may contain embedded typefaces. In accordance with Adobe's licensing policy, this file may be printed or viewed but
shall not be edited unless the typefaces which are embedded are licensed to and installed on the computer performing the editing. In
downloading this file, parties accept therein the responsibility of not infringing Adobe's licensing policy. The ISO Central Secretariat
accepts no liability in this area.
Adobe is a trademark of Adobe Systems Incorporated.
Details of the software products used to create this PDF file can be found in the General Info relative to the file; the PDF-creation
parameters were optimized for printing. Every care has been taken to ensure that the file is suitable for use by ISO member bodies. In
the unlikely event that a problem relating to it is found, please inform the Central Secretariat at the address given below.
COPYRIGHT PROTECTED DOCUMENT
© ISO 2010
All rights reserved. Unless otherwise specified, no part of this publication may be reproduced or utilized in any form or by any means,
electronic or mechanical, including photocopying and microfilm, without permission in writing from either ISO at the address below or
ISO's member body in the country of the requester.
ISO copyright office
Case postale 56 • CH-1211 Geneva 20
Tel. + 41 22 749 01 11
Fax + 41 22 749 09 47
E-mail copyright@iso.org
Web www.iso.org
Published in Switzerland
ii © ISO 2010 – All rights reserved
---------------------- Page: 2 ----------------------
ISO 28640:2010(E)
Contents Page
Foreword .iv
Introduction.v
1 Scope.1
2 Normative references.1
3 Terms and definitions .1
4 Symbols and mathematical binary operations.2
4.1 Symbols.2
4.2 Mathematical binary operations .2
5 Uniformly distributed pseudo-random numbers.3
5.1 General .3
5.2 M-sequence method definition .3
5.3 Pentanomial GFSR method .4
5.4 Combined Tausworthe method.4
5.5 Mersenne Twister method .5
6 Generation of random numbers from various distributions.6
6.1 Introduction.6
6.2 Uniform distribution .6
6.3 Standard beta distribution.7
6.4 Triangular distribution .8
6.5 General exponential distribution with location and scale parameters .8
6.6 Normal distribution .9
6.7 Gamma distribution.9
6.8 Weibull distribution .11
6.9 Lognormal distribution .11
6.10 Logistic distribution .11
6.11 Multivariate normal random variate generation .12
6.12 Binomial distribution.12
6.13 Poisson distribution.14
6.14 Discrete uniform distribution .14
Annex A (informative) Table of physical random numbers .16
A.1 Table of random numbers .16
A.2 Method of physical random number generation.17
Annex B (informative) Algorithm for pseudo-random number generation.18
B.1 Program code for the trinomial GFSR method.18
B.2 Program code for the pentanomial GFSR method.22
B.3 Program code for the combined Tausworthe method.28
B.4 Program code for the Mersenne Twister method .32
B.5 Linear congruential method .38
B.6 Reference examples.52
Bibliography.53
© ISO 2010 – All rights reserved iii
---------------------- Page: 3 ----------------------
ISO 28640:2010(E)
Foreword
ISO (the International Organization for Standardization) is a worldwide federation of national standards bodies
(ISO member bodies). The work of preparing International Standards is normally carried out through ISO
technical committees. Each member body interested in a subject for which a technical committee has been
established has the right to be represented on that committee. International organizations, governmental and
non-governmental, in liaison with ISO, also take part in the work. ISO collaborates closely with the
International Electrotechnical Commission (IEC) on all matters of electrotechnical standardization.
International Standards are drafted in accordance with the rules given in the ISO/IEC Directives, Part 2.
The main task of technical committees is to prepare International Standards. Draft International Standards
adopted by the technical committees are circulated to the member bodies for voting. Publication as an
International Standard requires approval by at least 75 % of the member bodies casting a vote.
Attention is drawn to the possibility that some of the elements of this document may be the subject of patent
rights. ISO shall not be held responsible for identifying any or all such patent rights.
ISO 28640 was prepared by Technical Committee ISO/TC 69, Applications of statistical methods.
This is the first edition.
iv © ISO 2010 – All rights reserved
---------------------- Page: 4 ----------------------
ISO 28640:2010(E)
Introduction
This International Standard specifies typical algorithms by which the users can regard the generated
numerical sequences as if they were real random variates.
Nowadays most statisticians, scientists and engineers have enough computer power at their disposal to carry
out large computer simulations, and it is important that these be based on sound pseudo-random generators.
This International Standard has been developed to help ensure that randomization, where needed, is carried
out correctly and efficiently.
Six uses of randomization can be identified in statistical standardization:
⎯ selection of a random sample;
⎯ analysis of sample data;
⎯ development of standards;
⎯ checking theoretical results;
⎯ demonstrating that a proposed procedure has the properties claimed of it;
⎯ resolving uncertainty in the statistical literature.
© ISO 2010 – All rights reserved v
---------------------- Page: 5 ----------------------
INTERNATIONAL STANDARD ISO 28640:2010(E)
Random variate generation methods
1 Scope
This International Standard specifies methods for generating uniform and non-uniform random variates for
Monte Carlo simulation purposes. Cryptographic random number generation methods are not included. This
International Standard is applicable, inter alia, by
⎯ researchers, industrial engineers or experts in operations management, who use statistical simulation,
⎯ statisticians who need randomization related to SQC methods, statistical design of experiments or sample
surveys,
⎯ applied mathematicians who plan complex optimization procedures that require the use of Monte Carlo
methods, and
⎯ software engineers who implement algorithms for random variate generation.
2 Normative references
The following referenced documents are indispensable for the application of this document. For dated
references, only the edition cited applies. For undated references, the latest edition of the referenced
document (including any amendments) applies.
ISO/IEC 2382-1, Information technology — Vocabulary — Part 1: Fundamental terms
ISO 3534-1, Statistics — Vocabulary and symbols — Part 1: General statistical terms and terms used in
probability
ISO 3534-2, Statistics — Vocabulary and symbols — Part 2: Applied statistics
3 Terms and definitions
For the purposes of this document, the terms and definitions given in ISO/IEC 2382-1, ISO 3534-1 and
ISO 3534-2 apply, except where redefined below.
3.1
random variate
random number
number as the realization of a specific random variable
NOTE 1 The term “random number” is often used for uniformly distributed random variate.
NOTE 2 Random numbers provided as a sequence are called a “random number sequence”.
© ISO 2010 – All rights reserved 1
---------------------- Page: 6 ----------------------
ISO 28640:2010(E)
3.2
pseudo-random number
random number (3.1) generated by an algorithm, that appears to be random
NOTE If there is no fear of misunderstanding, a pseudo-random number may simply be called a “random number”.
3.3
physical random number
random number (3.1) generated by a physical mechanism
3.4
binary random number sequence
random number (3.1) sequence consisting of zeros and ones
3.5
seed
initialization value required for pseudo-random number generation
4 Symbols and mathematical binary operations
4.1 Symbols
For the purposes of this document, the symbols given in the normative references as the latest versions of
ISO/IEC 2382-1, ISO 3534-1 and ISO 3534-2 apply, except where redefined below.
The symbols and abbreviations specifically used in this International Standard are as follows:
X integer type uniform random number
U standard uniform random number
Z normal random variate
n suffix of random number sequence
4.2 Mathematical binary operations
The mathematical binary operations specifically used in this International Standard are as follows:
mod(m; k) residue from dividing integer m by k
m ⊕ k bitwise exclusive logical disjunction of binary integers m and k
EXAMPLE 1 1 ⊕ 1 = 0
0 ⊕ 1 = 1
1 ⊕ 0 = 1
0 ⊕ 0 = 0
1010 ⊕ 1100 = 0110
2 © ISO 2010 – All rights reserved
---------------------- Page: 7 ----------------------
ISO 28640:2010(E)
m ∧ k bitwise logical conjunction of binary integers m and k
EXAMPLE 2 1 ∧ 1 = 1
0 ∧ 1 = 0
1 ∧ 0 = 0
0 ∧ 0 = 0
1010 ∧ 1100 = 1000
m := k replaces value m by k
m >> k k-bit right shift of binary integer m
m << k k-bit left shift of binary integer m
5 Uniformly distributed pseudo-random numbers
5.1 General
This clause provides algorithms for generating uniformly distributed pseudo-random numbers based on
M-sequence methods (see 5.2).
Annex A introduces the concept of physically generated random numbers for information.
Annex B includes C and full Basic codes for all the recommended algorithms for information. Although the
linear congruential method is not recommended for complex Monte Carlo simulations, it is also included in
Annex B for information.
5.2 M-sequence method definition
a) Let p be a natural number, and c , c , ., c be specified to be 0 or 1, and define the recurrence
1 2 p − 1
formula
x = c x + c x + . + c x + x (mod 2) (n = 1, 2, 3, .)
n + p p − 1 n + p − 1 p − 2 n + p − 2 1 n + 1 n
b) The least positive integer N such that x = x for all n is called the period of the sequence. This
n + N n
p
sequence is called an M-sequence in cases where its period is 2 − 1.
c) The polynomial
p p − 1
t + c t + . + c t + 1
p − 1 1
is called the characteristic polynomial of the above-mentioned recurrence formula.
NOTE 1 A necessary and sufficient condition for the above-mentioned recurrence formula to generate an M-sequence
is that at least one of the seeds x , x , ., x is not zero.
1 2 p
NOTE 2 The letter M of the M-sequence originates from the English word “maximum”, which means the largest. The
p
period of any sequence generated by the above recurrence formula cannot exceed 2 − 1. Therefore, if there is a series
p
that has a period of 2 − 1, it is the series that has the largest period.
NOTE 3 When this method is used, either one of the polynomials listed in Table 1 or another primitive polynomial listed
in the literature is chosen as the characteristic polynomial and its coefficients are used to define the recurrence formula
in a).
© ISO 2010 – All rights reserved 3
---------------------- Page: 8 ----------------------
ISO 28640:2010(E)
5.3 Pentanomial GFSR method
This method uses a characteristic polynomial of 5 terms, and it generates binary integer sequences of w bits
by the following recurrence formula. This algorithm is called the GFSR or “generalized feedback shift register”
random number generator.
X = X ⊕ X ⊕ X ⊕ X (n = 1, 2, 3, .)
n + p n + q1 n + q2 n + q3 n
The parameters are (p, q1, q2, q3, w) and X , ., X are initially given as seeds. Examples of parameters p, q ,
1 p 1
p
q , q giving the largest period 2 − 1 are indicated in Table 1.
2 3
Table 1 — Pentanomial characteristic polynomials
p q q q
1 2 3
89 20 40 69
107 31 57 82
127 22 63 83
521 86 197 447
607 167 307 461
1 279 339 630 988
2 203 585 1 197 1 656
2 281 577 1 109 1 709
3 217 809 1 621 2 381
4 253 1 093 2 254 3 297
4 423 1 171 2 273 3 299
9 689 2 799 5 463 7 712
NOTE q , q , q represent exponents of the non-zero terms of the
1 2 3
characteristic polynomial.
5.4 Combined Tausworthe method
Let x , x , x , … be an M-sequence generated by the recurrence relationship:
0 1 2
x = x + x (mod 2) (n = 0, 1, 2, …)
n + p n + q n
Using this M-sequence, a w-bit integer sequence called a simple Tausworthe sequence with parameters
(p, q, t) is obtained as follows:
X = x x …x (n = 0, 1, 2, …)
n nt nt + 1 nt + w − 1
where
p
t is a natural number which is coprime to the period 2 − 1 of the M-sequence;
w is the word length which does not exceed p.
p
The period of this sequence is also 2 − 1.
NOTE 1 Two integers are said to be coprime, or relatively prime, when they have no common divisors other than unity.
4 © ISO 2010 – All rights reserved
---------------------- Page: 9 ----------------------
ISO 28640:2010(E)
4
EXAMPLE If a primitive polynomial t + t + 1 is chosen, set p = 4, and q = 1 in the above recurrence relationship. If
the seeds (x , x , x , x ) = (1,1,1,1) are given to the recurrence, then the M-sequence obtained by the recurrence will be
0 1 2 3
4
1,1,1,1, 0,0,0,1, 0,0,1,1, 0,1,0,1, 1,1,1,0, … , and the period of the sequence is 2 − 1 = 15. Taking, for example, t = 4
which is coprime to 15, and w = 4, the simple Tausworthe sequence {X } with parameters (4, 1, 4) is obtained as follows:
n
X = x x x x = 1111 (= 15)
0 0 1 2 3
X = x x x x = 0001 (= 1)
1 4 5 6 7
X = x x x x = 0011 (= 3)
2 8 9 10 11
X = x x x x = 0101 (= 5)
3 12 13 14 0
X = x x x x = 1110 (= 14)
4 1 2 3 4
X = x x x x = 0010 (= 2)
5 5 6 7 8
.....
The simple Tausworthe sequence obtained in this way will be, in decimal notation, 15, 1, 3, 5, 14, 2, 6, 11, 12,
4
4, 13, 7, 8, 9, 10, 15, 1, 3, … , and its period is 2 − 1 = 15.
(j)
Suppose now that there is a multiple, say J, of simple Tausworthe sequences {X }, j = 1, 2, ., J with the
n
same word length w. The combined Tausworthe method is a technique that generates a sequence of
pseudo-random numbers {X } as the bitwise exclusive logical disjunction in the binary representation of these
n
J sequences.
(1) (2) (J)
X = X ⊕ X ⊕ … ⊕ X (n = 0, 1, 2, …)
n n n n
The parameters and the seeds of the combined Tausworthe sequence are combinations of the parameters
and the seeds of each simple Tausworthe sequence. If the periods of the J simple Tausworthe sequences are
coprime, then the period of the combined Tausworthe sequence is the product of the periods of the J
sequences.
NOTE 2 This method can generate sequences with good multidimensional equidistribution characteristics. The
algorithm taus88_31( ) given in Annex A generates a sequence of 31-bit integers by combining three simple Tausworthe
generators with parameters (p, q, t) = (31, 13, 12), (29, 2, 4), and (28, 3, 17), respectively. The period length of the
31 29 28 88
combined sequence is (2 − 1)(2 − 1)(2 − 1), i.e. about 2 . Many other combinations are suggested in References [7]
and [8] in the Bibliography.
5.5 Mersenne Twister method
Let X be a binary integer of w bits. Then, the Mersenne Twister method generates a sequence of binary
n
integer pseudo-random numbers of w bits according to the following recurrence formula with integer constants
p, q, r and a binary integer a of w bits.
f l (r)
X = X ⊕ (X |X ) A , (n = 1, 2, 3, .)
n + p n + q n n+1
f l (r) f l
where (X |X ) represents a binary integer that is obtained by a concatenation of X and X , the first
n n+1 n n + 1
w − r bits of X and the last r bits of X in this order. A is a w × w 0-1 matrix, which is determined by a, and
n n + 1
the product XA is given by the following formula.
X >> 1 (when the last bit of X = 0)
XA = (X >> 1) ⊕ a (when the last bit of X = 1)
Here, X is regarded as a w dimensional 0-1 vector.
© ISO 2010 – All rights reserved 5
---------------------- Page: 10 ----------------------
ISO 28640:2010(E)
pw−r
NOTE The necessary amount of memory for this computation is p words, the period becomes 2 − 1, and the
efficiency is better than that of the GFSR methods described previously. To improve the randomness of the first w − r bits,
the following series of conversions can be applied to X .
n
y := X
n
y := y ⊕ (y >> u)
y := y ⊕ [(y << s) ∧ b]
y := y ⊕ [(y << t) ∧ c]
y := y ⊕ (y >> l)
where b, c are constant bits masks to improve the randomness of the first w − r bits. The parameters of this
algorithm are (p, q, r, w, a, u, s, t, l, b, c). The seeds are X , ., X and the first w − r bits of X .
2 q + 1 1
The final value of y is the pseudo-random number.
6 Generation of random numbers from various distributions
6.1 Introduction
Methods of generating random numbers Y from various distributions by using uniform random numbers X of
integer type, are described below.
The distribution function is denoted by F(y). If it is a continuous distribution, its probability density function is
denoted by f(y), and if it is a discrete distribution, its probability mass function is denoted by p(y).
6.2 Uniform distribution
6.2.1 Standard uniform distribution
6.2.1.1 Probability density function
⎧1, 0uuy 1
fy() =
⎨
0, otherwise
⎩
6.2.1.2 Random variate generation method
If the maximum value of uniform random number X of integer type is m − 1, the following formula should be
used to generate standard uniform random numbers.
X
U =
m
EXAMPLE For any w-bit integer sequences generated by the method described in 5.2 through 5.5, m = 2.
NOTE 1 Because X takes on discrete values, the values of U are also discrete.
NOTE 2 The value of U never becomes 1,0. The value of U becomes 0,0 only when X = 0. In the case of M-sequence
random numbers, any generation method may cause this phenomenon.
NOTE 3 Random numbers from the standard uniform distribution are called standard uniform random numbers, and
are represented by U , U , … They are assumed to be independent of each other.
1 2
6 © ISO 2010 – All rights reserved
---------------------- Page: 11 ----------------------
ISO 28640:2010(E)
6.2.2 General uniform distribution
6.2.2.1 Probability density function
⎧1/ba, uuy a +b
fy() =
⎨
0, otherwise
⎩
where b > 0.
6.2.2.2 Random variate generation method
If the standard uniform random number U is generated by the method specified in 6.2.1.2, then the general
uniform random number should be generated by the following formula:
Yb=+U a
6.3 Standard beta distribution
6.3.1 Probability density function
d −1
⎧ c−1
yy1−
()
⎪
01uuy
,
fy() =
⎨
Β cd,
()
⎪
otherwise
0,
⎩
1
d −1
c−1
where Β cd,1=−x xdx is the beta function and the parameters c and d are greater than 0.
() ( )
∫
0
6.3.2 Random variate generation method by Jöhnk
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the standard beta random number Y should be generated by the following procedures.
1/cd1/ 1/ c
If YU=+U is less than or equal to 1, set YU= /Y ; otherwise, generate two standard uniform
12 1
random numbers again until the inequality is satisfied.
6.3.3 Random variate generation method by Cheng
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the standard beta random number Y should be generated by the following procedures.
[Set-up]
a) Let
⎧ < 1
min()cd, , if min(c,d)
⎪
q =
⎨
2(cd−+c d)
, otherwise
⎪
cd+− 2
⎩
[Generation]
b) Let
1 U
1
VW==,ecxp(V)
qU1−
1
© ISO 2010 – All rights reserved 7
---------------------- Page: 12 ----------------------
ISO 28640:2010(E)
c) If
⎛⎞cd+
2
()cd++ln (c+q)V−ln4Wln(UU)
⎜⎟ 12
dW+
⎝⎠
then employ
W
Y = ; and stop.
dW+
d) Generate U , U , and go to b).
1 2
U
1 ⎛⎞cd+ W
2
1
Vc=+()dln +(c+q)V−ln4Wln(UU)
⎜⎟ 12
q 1−+U dW dW+
⎝⎠
1
Jöhnk's method is recommended when max(c, d) u 1; otherwise, Cheng's method is recommended.
NOTE General beta random variates with the support [a, a + b] will be obtained by a linear transformation similar to
the one described in 6.2.2.2.
6.4 Triangular distribution
6.4.1 Probability density function
⎧
ba−−y
⎪
, ab−+uuy a b
fy() = 2
⎨
b
⎪
0, otherwise
⎩
where b > 0.
6.4.2 Random variate generation method
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the triangular random number Y should be generated by Ya=+b()U +U −1 .
12
6.5 General exponential distribution with location and scale parameters
6.5.1 Probability density function
1
⎧
exp{}−−()ya/b , yWa
⎪
fy() = b
⎨
⎪
⎩0, ya<
where a and b are the location and scale parameters of the exponential distribution, respectively.
6.5.2 Random variate generation method
If the standard uniform random number U is generated by the method specified in 6.2.1, then the general
exponential random number should be generated by
Y = − b ln(U) + a
8 © ISO 2010 – All rights reserved
---------------------- Page: 13 ----------------------
ISO 28640:2010(E)
6.6 Normal distribution
6.6.1 Probability density function
⎧⎫
11 2
fz=−exp z− µ ,−∞
() ()
⎨⎬
2
2πσ
⎩⎭2σ
where µ and σ are the mean and standard deviation of the normal distribution, respectively.
NOTE The symbol Z is used for a normal random variate.
6.6.2 The Box-Müller method
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then two independent normal random numbers Z , Z will be generated by the following procedures:
1 2
Z =+µσ −2ln 1−UUcos(2π )
()
11 2
Z 2=+µσ −2ln 1−UUsin(2π )
()
12
NOTE 1 Since U is not continuous, Z , Z cannot be normally distributed in a strict sense. For example, using this
1 1 2
procedure, the upper bound of the absolute value of the pseudo-standardized normal variates is
−1
32 31
M=−2ln()mm=2In ; thus, when m = 2 , M ≈ 6,660 4, and when m = 2 − 1, M ≈ 6,555 5. However, since the
−10
probability that absolute values of random variables from a true standard normal distribution exceed M is only about 10 ,
this will rarely be a problem in practice.
NOTE 2 When generating U , U , by a linear congruential method sequentially, U and U are not independent of each
1 2 1 2
other, so the tail of the distribution of the generated Z and Z can depart considerably from the true normal distribution.
1 2
6.7 Gamma distribution
6.7.1 Probability density function
⎧ 1
c−1
()ya−−/b exp(y−a)/b , yWa
{} { }
⎪
fy() = bcΓ()
⎨
⎪
se
⎩0, otherwi
where a, b, c are the location, scale and shape parameters of the distribution, respectively.
6.7.2 Random variate generation methods
6.7.2.1 General
Algorithms are given for three special cases depending on the shape parameter value c as follows.
6.7.2.2 Case of c = k (k: integer)
Using independent uniform random numbers U , U , . , U , the transformation formula
1 2 k
Ya=−b ln{(1−U )(1−U ).(1−U )}
12 k
should be used.
© ISO 2010 – All rights reserved 9
---------------------- Page: 14 ----------------------
ISO 28640:2010(E)
NOTE The chi-squared distribution with even degrees of freedom can be generated by this method when a = 0 and
b = 2.
6.7.2.3 Case of c = k + 1/2 (k: integer)
Using a standard normal random number Z and uniform random number U , U , . , U , the transformation
0 1 2 k
formula
2
⎡⎤
Ya=+bZ /2l−n (1−U )(1−U ).(1−U )
{}
012 k
⎣⎦
should be used. In the case where k = 0, the logarithm term disappears.
NOTE The chi-squared distribution with odd degrees of freedom can be generated by this method when a = 0 and
b = 2.
6.7.2.4 Approximate generation method when c > 1/3
3
a) Set rc=−1/3, , s= rt=r−rln()r, p=1/(3 s) and qr=−3 .
b) Generate standard normal random number Z.
c) If Z < q, then go to b).
32
d) Set Yp=+()Z s, 2V=Z/, and generate U.
2
e) If (Y − r) /Y − V u U, then employ Y := a + bY and stop.
f) Set W = Y − r ln(Y) − t − V.
g) If W u U, then employ Y := a + bY and stop.
h) If W > −In(1,0 − U), then go to b).
NOTE This method is based on the Wilson-Hilferty transformation of chi-square distributions to an approximate
standard normal distribution. The accuracy of this approximation depends on the parameter values of c, and the
dependency is rather complicated. A very rough idea is as follows: the absolute difference between the percentage point
of the approximate distribution and the exact distribution is always less than 0,2.
6.7.2.5 Exact generation method when c > 1/2, by Cheng
a) Set q = c – ln 4 and rc=+ 2c−1.
b) Generate standard uniform random numbers U and U .
1 2
U
1 2
c) Set Vc==ln(),W c exp(U),Z=U U and R = q + rV − W.
112
1−U
1
d) If R W 4,5Z − (1 + In 4,5) then employ Y = a + bW and stop.
e) If R W InZ then employ Y = a + bW and stop.
f) Go to b).
10 © ISO 2010 – All rights reserved
---------------------- Page: 15 ----------------------
ISO 28640:2010(E)
p
⎛⎞⎛ ⎞
UU
2
11
pc=1/2 −1, lq=−cn4 r=+c 2c−1 q+prln −c W4,5UU −(1+ln4,5)
⎜⎟⎜ ⎟
()12
11−−UU
⎝⎠11⎝ ⎠
p
⎛⎞
U
1
Ya=+bc
⎜⎟
1−U
⎝⎠1
6.8 Weibull distribution
6.8.1 Probability distribution function
c
⎧
⎧⎫
ya−
⎪⎪⎛⎞
⎪1e−−xp , yaW
⎪⎨ ⎬
⎜⎟
Fy() = b
⎝⎠
⎨
⎪⎪
⎩⎭
⎪
0, ya<
⎪
⎩
where a, b and c are the location, scale and shape parameters of the distribution, respectively.
6.8.2 Random variate generation method
If the standard uniform rando
...
SLOVENSKI STANDARD
SIST ISO 28640:2010
01-julij-2010
0HWRGHJHQHULUDQMDQDNOMXþQLKVSUHPHQOMLYN
Random variate generation methods
Méthodes de génération de nombres pseudo-aléatoires
Ta slovenski standard je istoveten z: ISO 28640:2010
ICS:
03.120.30 8SRUDEDVWDWLVWLþQLKPHWRG Application of statistical
methods
SIST ISO 28640:2010 en
2003-01.Slovenski inštitut za standardizacijo. Razmnoževanje celote ali delov tega standarda ni dovoljeno.
---------------------- Page: 1 ----------------------
SIST ISO 28640:2010
---------------------- Page: 2 ----------------------
SIST ISO 28640:2010
INTERNATIONAL ISO
STANDARD 28640
First edition
2010-03-15
Random variate generation methods
Méthodes de génération de nombres pseudo-aléatoires
Reference number
ISO 28640:2010(E)
©
ISO 2010
---------------------- Page: 3 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
PDF disclaimer
This PDF file may contain embedded typefaces. In accordance with Adobe's licensing policy, this file may be printed or viewed but
shall not be edited unless the typefaces which are embedded are licensed to and installed on the computer performing the editing. In
downloading this file, parties accept therein the responsibility of not infringing Adobe's licensing policy. The ISO Central Secretariat
accepts no liability in this area.
Adobe is a trademark of Adobe Systems Incorporated.
Details of the software products used to create this PDF file can be found in the General Info relative to the file; the PDF-creation
parameters were optimized for printing. Every care has been taken to ensure that the file is suitable for use by ISO member bodies. In
the unlikely event that a problem relating to it is found, please inform the Central Secretariat at the address given below.
COPYRIGHT PROTECTED DOCUMENT
© ISO 2010
All rights reserved. Unless otherwise specified, no part of this publication may be reproduced or utilized in any form or by any means,
electronic or mechanical, including photocopying and microfilm, without permission in writing from either ISO at the address below or
ISO's member body in the country of the requester.
ISO copyright office
Case postale 56 • CH-1211 Geneva 20
Tel. + 41 22 749 01 11
Fax + 41 22 749 09 47
E-mail copyright@iso.org
Web www.iso.org
Published in Switzerland
ii © ISO 2010 – All rights reserved
---------------------- Page: 4 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
Contents Page
Foreword .iv
Introduction.v
1 Scope.1
2 Normative references.1
3 Terms and definitions .1
4 Symbols and mathematical binary operations.2
4.1 Symbols.2
4.2 Mathematical binary operations .2
5 Uniformly distributed pseudo-random numbers.3
5.1 General .3
5.2 M-sequence method definition .3
5.3 Pentanomial GFSR method .4
5.4 Combined Tausworthe method.4
5.5 Mersenne Twister method .5
6 Generation of random numbers from various distributions.6
6.1 Introduction.6
6.2 Uniform distribution .6
6.3 Standard beta distribution.7
6.4 Triangular distribution .8
6.5 General exponential distribution with location and scale parameters .8
6.6 Normal distribution .9
6.7 Gamma distribution.9
6.8 Weibull distribution .11
6.9 Lognormal distribution .11
6.10 Logistic distribution .11
6.11 Multivariate normal random variate generation .12
6.12 Binomial distribution.12
6.13 Poisson distribution.14
6.14 Discrete uniform distribution .14
Annex A (informative) Table of physical random numbers .16
A.1 Table of random numbers .16
A.2 Method of physical random number generation.17
Annex B (informative) Algorithm for pseudo-random number generation.18
B.1 Program code for the trinomial GFSR method.18
B.2 Program code for the pentanomial GFSR method.22
B.3 Program code for the combined Tausworthe method.28
B.4 Program code for the Mersenne Twister method .32
B.5 Linear congruential method .38
B.6 Reference examples.52
Bibliography.53
© ISO 2010 – All rights reserved iii
---------------------- Page: 5 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
Foreword
ISO (the International Organization for Standardization) is a worldwide federation of national standards bodies
(ISO member bodies). The work of preparing International Standards is normally carried out through ISO
technical committees. Each member body interested in a subject for which a technical committee has been
established has the right to be represented on that committee. International organizations, governmental and
non-governmental, in liaison with ISO, also take part in the work. ISO collaborates closely with the
International Electrotechnical Commission (IEC) on all matters of electrotechnical standardization.
International Standards are drafted in accordance with the rules given in the ISO/IEC Directives, Part 2.
The main task of technical committees is to prepare International Standards. Draft International Standards
adopted by the technical committees are circulated to the member bodies for voting. Publication as an
International Standard requires approval by at least 75 % of the member bodies casting a vote.
Attention is drawn to the possibility that some of the elements of this document may be the subject of patent
rights. ISO shall not be held responsible for identifying any or all such patent rights.
ISO 28640 was prepared by Technical Committee ISO/TC 69, Applications of statistical methods.
This is the first edition.
iv © ISO 2010 – All rights reserved
---------------------- Page: 6 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
Introduction
This International Standard specifies typical algorithms by which the users can regard the generated
numerical sequences as if they were real random variates.
Nowadays most statisticians, scientists and engineers have enough computer power at their disposal to carry
out large computer simulations, and it is important that these be based on sound pseudo-random generators.
This International Standard has been developed to help ensure that randomization, where needed, is carried
out correctly and efficiently.
Six uses of randomization can be identified in statistical standardization:
⎯ selection of a random sample;
⎯ analysis of sample data;
⎯ development of standards;
⎯ checking theoretical results;
⎯ demonstrating that a proposed procedure has the properties claimed of it;
⎯ resolving uncertainty in the statistical literature.
© ISO 2010 – All rights reserved v
---------------------- Page: 7 ----------------------
SIST ISO 28640:2010
---------------------- Page: 8 ----------------------
SIST ISO 28640:2010
INTERNATIONAL STANDARD ISO 28640:2010(E)
Random variate generation methods
1 Scope
This International Standard specifies methods for generating uniform and non-uniform random variates for
Monte Carlo simulation purposes. Cryptographic random number generation methods are not included. This
International Standard is applicable, inter alia, by
⎯ researchers, industrial engineers or experts in operations management, who use statistical simulation,
⎯ statisticians who need randomization related to SQC methods, statistical design of experiments or sample
surveys,
⎯ applied mathematicians who plan complex optimization procedures that require the use of Monte Carlo
methods, and
⎯ software engineers who implement algorithms for random variate generation.
2 Normative references
The following referenced documents are indispensable for the application of this document. For dated
references, only the edition cited applies. For undated references, the latest edition of the referenced
document (including any amendments) applies.
ISO/IEC 2382-1, Information technology — Vocabulary — Part 1: Fundamental terms
ISO 3534-1, Statistics — Vocabulary and symbols — Part 1: General statistical terms and terms used in
probability
ISO 3534-2, Statistics — Vocabulary and symbols — Part 2: Applied statistics
3 Terms and definitions
For the purposes of this document, the terms and definitions given in ISO/IEC 2382-1, ISO 3534-1 and
ISO 3534-2 apply, except where redefined below.
3.1
random variate
random number
number as the realization of a specific random variable
NOTE 1 The term “random number” is often used for uniformly distributed random variate.
NOTE 2 Random numbers provided as a sequence are called a “random number sequence”.
© ISO 2010 – All rights reserved 1
---------------------- Page: 9 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
3.2
pseudo-random number
random number (3.1) generated by an algorithm, that appears to be random
NOTE If there is no fear of misunderstanding, a pseudo-random number may simply be called a “random number”.
3.3
physical random number
random number (3.1) generated by a physical mechanism
3.4
binary random number sequence
random number (3.1) sequence consisting of zeros and ones
3.5
seed
initialization value required for pseudo-random number generation
4 Symbols and mathematical binary operations
4.1 Symbols
For the purposes of this document, the symbols given in the normative references as the latest versions of
ISO/IEC 2382-1, ISO 3534-1 and ISO 3534-2 apply, except where redefined below.
The symbols and abbreviations specifically used in this International Standard are as follows:
X integer type uniform random number
U standard uniform random number
Z normal random variate
n suffix of random number sequence
4.2 Mathematical binary operations
The mathematical binary operations specifically used in this International Standard are as follows:
mod(m; k) residue from dividing integer m by k
m ⊕ k bitwise exclusive logical disjunction of binary integers m and k
EXAMPLE 1 1 ⊕ 1 = 0
0 ⊕ 1 = 1
1 ⊕ 0 = 1
0 ⊕ 0 = 0
1010 ⊕ 1100 = 0110
2 © ISO 2010 – All rights reserved
---------------------- Page: 10 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
m ∧ k bitwise logical conjunction of binary integers m and k
EXAMPLE 2 1 ∧ 1 = 1
0 ∧ 1 = 0
1 ∧ 0 = 0
0 ∧ 0 = 0
1010 ∧ 1100 = 1000
m := k replaces value m by k
m >> k k-bit right shift of binary integer m
m << k k-bit left shift of binary integer m
5 Uniformly distributed pseudo-random numbers
5.1 General
This clause provides algorithms for generating uniformly distributed pseudo-random numbers based on
M-sequence methods (see 5.2).
Annex A introduces the concept of physically generated random numbers for information.
Annex B includes C and full Basic codes for all the recommended algorithms for information. Although the
linear congruential method is not recommended for complex Monte Carlo simulations, it is also included in
Annex B for information.
5.2 M-sequence method definition
a) Let p be a natural number, and c , c , ., c be specified to be 0 or 1, and define the recurrence
1 2 p − 1
formula
x = c x + c x + . + c x + x (mod 2) (n = 1, 2, 3, .)
n + p p − 1 n + p − 1 p − 2 n + p − 2 1 n + 1 n
b) The least positive integer N such that x = x for all n is called the period of the sequence. This
n + N n
p
sequence is called an M-sequence in cases where its period is 2 − 1.
c) The polynomial
p p − 1
t + c t + . + c t + 1
p − 1 1
is called the characteristic polynomial of the above-mentioned recurrence formula.
NOTE 1 A necessary and sufficient condition for the above-mentioned recurrence formula to generate an M-sequence
is that at least one of the seeds x , x , ., x is not zero.
1 2 p
NOTE 2 The letter M of the M-sequence originates from the English word “maximum”, which means the largest. The
p
period of any sequence generated by the above recurrence formula cannot exceed 2 − 1. Therefore, if there is a series
p
that has a period of 2 − 1, it is the series that has the largest period.
NOTE 3 When this method is used, either one of the polynomials listed in Table 1 or another primitive polynomial listed
in the literature is chosen as the characteristic polynomial and its coefficients are used to define the recurrence formula
in a).
© ISO 2010 – All rights reserved 3
---------------------- Page: 11 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
5.3 Pentanomial GFSR method
This method uses a characteristic polynomial of 5 terms, and it generates binary integer sequences of w bits
by the following recurrence formula. This algorithm is called the GFSR or “generalized feedback shift register”
random number generator.
X = X ⊕ X ⊕ X ⊕ X (n = 1, 2, 3, .)
n + p n + q1 n + q2 n + q3 n
The parameters are (p, q1, q2, q3, w) and X , ., X are initially given as seeds. Examples of parameters p, q ,
1 p 1
p
q , q giving the largest period 2 − 1 are indicated in Table 1.
2 3
Table 1 — Pentanomial characteristic polynomials
p q q q
1 2 3
89 20 40 69
107 31 57 82
127 22 63 83
521 86 197 447
607 167 307 461
1 279 339 630 988
2 203 585 1 197 1 656
2 281 577 1 109 1 709
3 217 809 1 621 2 381
4 253 1 093 2 254 3 297
4 423 1 171 2 273 3 299
9 689 2 799 5 463 7 712
NOTE q , q , q represent exponents of the non-zero terms of the
1 2 3
characteristic polynomial.
5.4 Combined Tausworthe method
Let x , x , x , … be an M-sequence generated by the recurrence relationship:
0 1 2
x = x + x (mod 2) (n = 0, 1, 2, …)
n + p n + q n
Using this M-sequence, a w-bit integer sequence called a simple Tausworthe sequence with parameters
(p, q, t) is obtained as follows:
X = x x …x (n = 0, 1, 2, …)
n nt nt + 1 nt + w − 1
where
p
t is a natural number which is coprime to the period 2 − 1 of the M-sequence;
w is the word length which does not exceed p.
p
The period of this sequence is also 2 − 1.
NOTE 1 Two integers are said to be coprime, or relatively prime, when they have no common divisors other than unity.
4 © ISO 2010 – All rights reserved
---------------------- Page: 12 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
4
EXAMPLE If a primitive polynomial t + t + 1 is chosen, set p = 4, and q = 1 in the above recurrence relationship. If
the seeds (x , x , x , x ) = (1,1,1,1) are given to the recurrence, then the M-sequence obtained by the recurrence will be
0 1 2 3
4
1,1,1,1, 0,0,0,1, 0,0,1,1, 0,1,0,1, 1,1,1,0, … , and the period of the sequence is 2 − 1 = 15. Taking, for example, t = 4
which is coprime to 15, and w = 4, the simple Tausworthe sequence {X } with parameters (4, 1, 4) is obtained as follows:
n
X = x x x x = 1111 (= 15)
0 0 1 2 3
X = x x x x = 0001 (= 1)
1 4 5 6 7
X = x x x x = 0011 (= 3)
2 8 9 10 11
X = x x x x = 0101 (= 5)
3 12 13 14 0
X = x x x x = 1110 (= 14)
4 1 2 3 4
X = x x x x = 0010 (= 2)
5 5 6 7 8
.....
The simple Tausworthe sequence obtained in this way will be, in decimal notation, 15, 1, 3, 5, 14, 2, 6, 11, 12,
4
4, 13, 7, 8, 9, 10, 15, 1, 3, … , and its period is 2 − 1 = 15.
(j)
Suppose now that there is a multiple, say J, of simple Tausworthe sequences {X }, j = 1, 2, ., J with the
n
same word length w. The combined Tausworthe method is a technique that generates a sequence of
pseudo-random numbers {X } as the bitwise exclusive logical disjunction in the binary representation of these
n
J sequences.
(1) (2) (J)
X = X ⊕ X ⊕ … ⊕ X (n = 0, 1, 2, …)
n n n n
The parameters and the seeds of the combined Tausworthe sequence are combinations of the parameters
and the seeds of each simple Tausworthe sequence. If the periods of the J simple Tausworthe sequences are
coprime, then the period of the combined Tausworthe sequence is the product of the periods of the J
sequences.
NOTE 2 This method can generate sequences with good multidimensional equidistribution characteristics. The
algorithm taus88_31( ) given in Annex A generates a sequence of 31-bit integers by combining three simple Tausworthe
generators with parameters (p, q, t) = (31, 13, 12), (29, 2, 4), and (28, 3, 17), respectively. The period length of the
31 29 28 88
combined sequence is (2 − 1)(2 − 1)(2 − 1), i.e. about 2 . Many other combinations are suggested in References [7]
and [8] in the Bibliography.
5.5 Mersenne Twister method
Let X be a binary integer of w bits. Then, the Mersenne Twister method generates a sequence of binary
n
integer pseudo-random numbers of w bits according to the following recurrence formula with integer constants
p, q, r and a binary integer a of w bits.
f l (r)
X = X ⊕ (X |X ) A , (n = 1, 2, 3, .)
n + p n + q n n+1
f l (r) f l
where (X |X ) represents a binary integer that is obtained by a concatenation of X and X , the first
n n+1 n n + 1
w − r bits of X and the last r bits of X in this order. A is a w × w 0-1 matrix, which is determined by a, and
n n + 1
the product XA is given by the following formula.
X >> 1 (when the last bit of X = 0)
XA = (X >> 1) ⊕ a (when the last bit of X = 1)
Here, X is regarded as a w dimensional 0-1 vector.
© ISO 2010 – All rights reserved 5
---------------------- Page: 13 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
pw−r
NOTE The necessary amount of memory for this computation is p words, the period becomes 2 − 1, and the
efficiency is better than that of the GFSR methods described previously. To improve the randomness of the first w − r bits,
the following series of conversions can be applied to X .
n
y := X
n
y := y ⊕ (y >> u)
y := y ⊕ [(y << s) ∧ b]
y := y ⊕ [(y << t) ∧ c]
y := y ⊕ (y >> l)
where b, c are constant bits masks to improve the randomness of the first w − r bits. The parameters of this
algorithm are (p, q, r, w, a, u, s, t, l, b, c). The seeds are X , ., X and the first w − r bits of X .
2 q + 1 1
The final value of y is the pseudo-random number.
6 Generation of random numbers from various distributions
6.1 Introduction
Methods of generating random numbers Y from various distributions by using uniform random numbers X of
integer type, are described below.
The distribution function is denoted by F(y). If it is a continuous distribution, its probability density function is
denoted by f(y), and if it is a discrete distribution, its probability mass function is denoted by p(y).
6.2 Uniform distribution
6.2.1 Standard uniform distribution
6.2.1.1 Probability density function
⎧1, 0uuy 1
fy() =
⎨
0, otherwise
⎩
6.2.1.2 Random variate generation method
If the maximum value of uniform random number X of integer type is m − 1, the following formula should be
used to generate standard uniform random numbers.
X
U =
m
EXAMPLE For any w-bit integer sequences generated by the method described in 5.2 through 5.5, m = 2.
NOTE 1 Because X takes on discrete values, the values of U are also discrete.
NOTE 2 The value of U never becomes 1,0. The value of U becomes 0,0 only when X = 0. In the case of M-sequence
random numbers, any generation method may cause this phenomenon.
NOTE 3 Random numbers from the standard uniform distribution are called standard uniform random numbers, and
are represented by U , U , … They are assumed to be independent of each other.
1 2
6 © ISO 2010 – All rights reserved
---------------------- Page: 14 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
6.2.2 General uniform distribution
6.2.2.1 Probability density function
⎧1/ba, uuy a +b
fy() =
⎨
0, otherwise
⎩
where b > 0.
6.2.2.2 Random variate generation method
If the standard uniform random number U is generated by the method specified in 6.2.1.2, then the general
uniform random number should be generated by the following formula:
Yb=+U a
6.3 Standard beta distribution
6.3.1 Probability density function
d −1
⎧ c−1
yy1−
()
⎪
01uuy
,
fy() =
⎨
Β cd,
()
⎪
otherwise
0,
⎩
1
d −1
c−1
where Β cd,1=−x xdx is the beta function and the parameters c and d are greater than 0.
() ( )
∫
0
6.3.2 Random variate generation method by Jöhnk
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the standard beta random number Y should be generated by the following procedures.
1/cd1/ 1/ c
If YU=+U is less than or equal to 1, set YU= /Y ; otherwise, generate two standard uniform
12 1
random numbers again until the inequality is satisfied.
6.3.3 Random variate generation method by Cheng
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the standard beta random number Y should be generated by the following procedures.
[Set-up]
a) Let
⎧ < 1
min()cd, , if min(c,d)
⎪
q =
⎨
2(cd−+c d)
, otherwise
⎪
cd+− 2
⎩
[Generation]
b) Let
1 U
1
VW==,ecxp(V)
qU1−
1
© ISO 2010 – All rights reserved 7
---------------------- Page: 15 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
c) If
⎛⎞cd+
2
()cd++ln (c+q)V−ln4Wln(UU)
⎜⎟ 12
dW+
⎝⎠
then employ
W
Y = ; and stop.
dW+
d) Generate U , U , and go to b).
1 2
U
1 ⎛⎞cd+ W
2
1
Vc=+()dln +(c+q)V−ln4Wln(UU)
⎜⎟ 12
q 1−+U dW dW+
⎝⎠
1
Jöhnk's method is recommended when max(c, d) u 1; otherwise, Cheng's method is recommended.
NOTE General beta random variates with the support [a, a + b] will be obtained by a linear transformation similar to
the one described in 6.2.2.2.
6.4 Triangular distribution
6.4.1 Probability density function
⎧
ba−−y
⎪
, ab−+uuy a b
fy() = 2
⎨
b
⎪
0, otherwise
⎩
where b > 0.
6.4.2 Random variate generation method
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the triangular random number Y should be generated by Ya=+b()U +U −1 .
12
6.5 General exponential distribution with location and scale parameters
6.5.1 Probability density function
1
⎧
exp{}−−()ya/b , yWa
⎪
fy() = b
⎨
⎪
⎩0, ya<
where a and b are the location and scale parameters of the exponential distribution, respectively.
6.5.2 Random variate generation method
If the standard uniform random number U is generated by the method specified in 6.2.1, then the general
exponential random number should be generated by
Y = − b ln(U) + a
8 © ISO 2010 – All rights reserved
---------------------- Page: 16 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
6.6 Normal distribution
6.6.1 Probability density function
⎧⎫
11 2
fz=−exp z− µ ,−∞
() ()
⎨⎬
2
2πσ
⎩⎭2σ
where µ and σ are the mean and standard deviation of the normal distribution, respectively.
NOTE The symbol Z is used for a normal random variate.
6.6.2 The Box-Müller method
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then two independent normal random numbers Z , Z will be generated by the following procedures:
1 2
Z =+µσ −2ln 1−UUcos(2π )
()
11 2
Z 2=+µσ −2ln 1−UUsin(2π )
()
12
NOTE 1 Since U is not continuous, Z , Z cannot be normally distributed in a strict sense. For example, using this
1 1 2
procedure, the upper bound of the absolute value of the pseudo-standardized normal variates is
−1
32 31
M=−2ln()mm=2In ; thus, when m = 2 , M ≈ 6,660 4, and when m = 2 − 1, M ≈ 6,555 5. However, since the
−10
probability that absolute values of random variables from a true standard normal distribution exceed M is only about 10 ,
this will rarely be a problem in practice.
NOTE 2 When generating U , U , by a linear congruential method sequentially, U and U are not independent of each
1 2 1 2
other, so the tail of the distribution of the generated Z and Z can depart considerably from the true normal distribution.
1 2
6.7 Gamma distribution
6.7.1 Probability density function
⎧ 1
c−1
()ya−−/b exp(y−a)/b , yWa
{} { }
⎪
fy() = bcΓ()
⎨
⎪
se
⎩0, otherwi
where a, b, c are the location, scale and shape parameters of the distribution, respectively.
6.7.2 Random variate generation methods
6.7.2.1 General
Algorithms are given for three special cases depending on the shape parameter value c as follows.
6.7.2.2 Case of c = k (k: integer)
Using independent uniform random numbers U , U , . , U , the transformation formula
1 2 k
Ya=−b ln{(1−U )(1−U ).(1−U )}
12 k
should be used.
© ISO 2010 – All rights reserved 9
---------------------- Page: 17 ----------------------
SIST ISO 28640:2010
ISO 28640:2010(E)
NOTE The chi-squared distribution with even degrees of freedom can be generated by this method when a = 0 and
b = 2.
6.7.2.3 Case of c = k + 1/2 (k: integer)
Using a standard normal random number Z and uniform random number U , U , . , U , the transformation
0 1 2 k
formula
2
⎡⎤
Ya=+bZ /2l−n (1−U )(1−U ).(1−U )
{}
012 k
⎣⎦
should be used. In the case where k = 0, the logarithm term disappears.
NOTE The chi-squared distribution with odd degrees of freedom can be generated by this method when a = 0 and
b = 2.
6.7.2.4 Approximate generation method when c > 1/3
3
a) Set rc=−1/3, , s= rt=r−rln()r, p=1/(3 s) and qr=−3 .
b) Generate standard normal random number Z.
c) If Z < q, then go to b).
32
d) Set Yp=+()Z s, 2V=Z/, and generate U.
2
e) If (Y − r) /Y − V u U, then employ Y := a + bY and stop.
f) Set W = Y − r ln(Y) − t − V.
g) If W u U, then employ Y := a + bY and stop.
h) If W > −In(1,0 − U), then go to b).
NOTE This method is based on the Wilson-Hilferty transformation of chi-square distributions to an approximate
standard normal distribution. The accuracy of this approximation depends on the parameter values of c, and the
dependency is rather complicated. A very rough idea is as follows: the absolute difference between the percentage point
of the approxi
...
Questions, Comments and Discussion
Ask us and Technical Secretary will try to provide an answer. You can facilitate discussion about the standard in here.